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1. Introduction 

The quest for finding optimum solutions to engineering problems has existed for a 
long time. In modern times, the development of optimization as a branch of applied 
mathematics is regarded to have originated in the works of Newton, Bernoulli and Euler. 
Venkayya has presented a historical perspective on optimization in [1]. The term 
'optimization' is defined by Ashley [2] as a procedure "...which attempts to choose the 
variables in a design process so as formally to achieve the best value of some 
performance index while not violating any of the associated conditions or constraints". 
Ashley presented an extensive review of practical applications of optimization in the 
aeronautical field till about 1980 [2]. It was noted that there existed an enormous amount 
of published literature in the field of optimization, but its practical applications in 
industry were very limited. Over the past 15 years, though, optimization has been widely 
applied to address practical problems in aerospace design [3-5]. 

The design of high performance aerospace systems is a complex task. It involves the 
integration of several disciplines such as aerodynamics, structural analysis, dynamics, and 
aeroelasticity. The problem involves multiple objectives and constraints pertaining to the 
design criteria associated with each of these disciplines. Many important trade-offs exist 
between the parameters involved which are used to define the different disciplines. 
Therefore, the development of multidisciplinary design optimization (MDO) techniques, 
in which different disciplines and design parameters are coupled into a closed loop 
numerical procedure, seems appropriate to address such a complex problem. The 
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importance of MDO in successful design of aerospace systems has been long recognized. 
Recent developments in this field have been surveyed by Sobieszczanski-Sobieski and 
Haftka [6]. 


1.1 Multidisciplinary Design Optimization of Smart Composite Wings 

The use of multidisciplinary optimization techniques in aircraft design has been 
increasing over the past decades. A recent survey of MDO applications in preliminary 
aircraft design has been presented by Kroo [7] . The analysis methods used in MDO 
range from simple analytical or empirical expressions to complex finite element models 
[8-11], The validity of the designs obtained using MDO procedures depends strongly 
upon the accuracy of the analytical methods used. Therefore, it is essential to integrate 
accurate and efficient analysis techniques to obtain meaningful optimum designs. 

Due to the high stiffness-to-weight ratio and directional stiffness and strength 
properties offered by advanced composites, they are increasingly being used in aircraft 
wing design. The deformation of the wing under load is directly related to the stacking 
sequence of the orthotropic composite laminae. Hence, for aircraft wings made out of 
composite materials, aeroelastic tailoring presents an opportunity to enhance structural, 
aerodynamic and control performance by utilizing their unique stiffness and strength 
properties. A comprehensive review of the aeroelastic tailoring technology was 
presented by Shirk et al. [12]. It was noted that aeroelastic tailoring has matured as a 
result of developments in composite material analysis and optimization techniques. Since 
aircraft wing design is a multidisciplinary problem involving the coupling of various 
disciplines such as aerodynamics, composite structural analysis, dynamics and 
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aeroelasticity, the use of MDO techniques seems appropriate. Key elements associated 
with the use of MDO in wing design are the development of sophisticated optimization 
techniques and accurate analysis methods. 


1.1.1 Optimization Techniques 

The design of composite wings to achieve high flutter and divergence speeds while 
maintaining a low structural weight and stresses within allowable limits, is a truly 
multidisciplinary problem. A number of optimization procedures have been developed 
for minimum weight design of wing structures with aeroelastic and other constraints. 
Most notable among these is the Automated Structural Optimization System (ASTROS) 
[13]. The procedure uses the finite element method for structural analysis and gradient 
based techniques for optimization. 

In general, an optimization problem can be associated with several objective 
functions, constraints and design variables. However, in many existing procedures the 
problem is formulated with a single objective function subject to several constraints [13]. 
Such procedures do not allow simultaneous minimization or maximization of more than 
one objective. A commonly used technique for addressing multiobjective problems is to 
combine individual objective functions in a linear fashion using weight factors [14]. 
Such methods are judgmental in nature as the answer depends upon weight factors which 
are often hard to justify. Also, the procedures do not satisfy the Kuhn-Tucker conditions 
of optimality [14]. To address this issue in a mathematically rigorous way, several 
formal multiobjective techniques have been developed by Chattopadhyay and McCarthy 
[15]. In the first method, called the modified global criterion approach, the individual 
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objective functions are combined into a single composite function using the optimized 
values for each objective function. The optimized values are obtained through single 
objective minimization or maximization prior to multiobjective formulation. The 
drawback of the method is that several single objective optimizations have to be 
performed to obtain the optimized values. Therefore, this method is not very efficient 
computationally. Another method, called 'Min (Minimum sum beta) [15], uses 
pseudo-design variables that represent deviations of the individual objective functions 
from respective target values. The objective function is then defined as a linear 
combination of these pseudo-design variables. An appealing feature of this approach is 
the fact that the objective function is a purely linear function. However, the prescription 
of the target values is rather arbitrary and hence prone to error. 

The multiobjective optimization formulation used in the present work is based on the 
Kreisselmeier-Steinhauser (K-S) function technique [16]. The K-S function approach can 
be applied to problems with multiple objectives and inequality constraints [17]. In this 
approach, each of the original objective functions is transformed into reduced objective 
functions using the value of the original objective functions calculated at the beginning of 
the cycle. The reduced objective functions are analogous to constraints. A new 
constraint vector is defined combining the original constraints and the reduced objective 
functions. Using the K-S function, a new objective function is now obtained which 
represents the original objective functions and constraints. Thus the K-S function 
technique efficiently integrates the objective functions and constraints into a single 
envelop function. The design variable vector in this formulation remains unchanged. 
The resulting unconstrained problem can be solved using any of the standard techniques 
for nonlinear optimization [14]. The search direction vector is obtained using the 
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Broyden-Fletcher-Goldberg-Shanno (BFGS) algorithm. With the search direction 
determined, a one-dimensional search for the minimum of the K-S function is performed 
using the three point quadratic approximation for step length [14], The gradients for 
objective functions and constraints are evaluated using the finite difference method. 

Since during an optimization process, several computations of the objective functions 
and constraints are necessary, it is computationally expensive to use exact analysis all the 
time. Therefore, approximation techniques are commonly used to reduce the analysis 
effort. Several different approximation techniques have been used in the literature, the 
most notable being the linear Taylor series and the reciprocal Taylor series 
approximations [14]. In this research, the two-point approximation technique [18], which 
has been found to be well suited for nonlinear optimization problems, is used. This 
technique uses the gradient information from the previous and current design cycles to 
construct the approximate function value. In the limiting cases the expansion reduces to 
the first order Taylor series or the reciprocal approximation form. 

In composite wings, the ply stacking sequence has a very strong influence on the 
design objectives and constraints. Often, ply angles are treated as continuous design 
variables, such as in ASTROS [13], and the resulting solution is replaced by the nearest 
integral value which can lead to suboptimal design. In practice, ply orientations (and ply 
thickness) are selected from a range of practical discrete values. Hence ply angles are 
best treated as discrete design variables. However, parameters such as wing span, chord 
and other dimensions represent continuous variables. Therefore, the composite wing 
optimization problem involves both continuous and discrete design variables. Recently, a 
hybrid optimization technique has been developed by Seeley and Chattopadhyay to 
simultaneously include continuous and discrete design variables in the optimization 
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problem [19], Since the hybrid optimization uses a combinatorial search technique [14] 
for discrete design variables, the procedure can be expensive computationally. Therefore, 
this approach has been used for only some of the optimization problems addressed in the 
present research. 

1.1.2 Composite Structural Analysis 

The structural analysis of aircraft wings can be performed either through a detailed 
investigation of the wing sections comprising skins, spars, ribs etc. or through the use of 
reduced structural models. The detailed analysis based on a full three-dimensional finite 
element solution [13] is computationally very expensive and requires a large modeling 
time. Hence, such techniques can be impractical in design optimization or trade-off 
studies during the conceptual design phase. Therefore, procedures based on reduced 
structural models, such as 'box beam' and 'equivalent plate', are frequently used during 
conceptual wing design. 

Among the aeroelastic analysis and optimization procedures based on reduced 
structural model, TSO [20] and ELAPS [21-22] have been widely used. These 
procedures use an equivalent plate model for structural analysis. The wing box geometry 
in TSO is limited to trapezoidal planforms, whereas ELAPS can analyze cranked wings 
through multiple trapezoidal segments. The depth of the structural box, which consists of 
cover skins and rib, stiffener and spar caps, can be varied over the planform using this 
procedure. Another procedure, named LS-CLASS [23], uses a structural model similar to 
ELAPS and includes analytical sensitivity derivatives for efficient aeroservoelastic 
optimization. However, all of these techniques are based on the Classical Laminate 
Theory (CLT) [24] which assumes that normals to the midplane before deformation 
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remain straight and normal to the plane after deformation. Thus CLT assumes 
deformation due entirely to bending and inplane stretching and neglects transverse shear 
stresses. Experimental results indicate that CLT underpredicts deflections and 
overpredicts natural frequencies [24]. Extensive studies comparing the results from CLT 
based equivalent plate models to detailed finite element models have shown serious 
limitations of the former approach [25-26]. The largest errors are reported in the torsional 
behavior. The representation of the transverse shear is identified as the principal reason 
for the differences. An equivalent plate model based on the First-order Shear 
Deformation Theory (FSDT) [24] yields a better correlation with finite element solution 
[25-26]. Thus, the inclusion of transverse shear in equivalent plate models is very 
important. The FSDT assumes constant transverse shear strain through the laminate 
thickness. This theory requires shear correction factors which are difficult to obtain 
because they depend on the lamina properties and the lamination scheme. Therefore, a 
more accurate description of the transverse shear stresses is necessary. 

The Higher-Order Shear Deformation Theory (HSDT) is capable of accurately and 
efficiently predicting the transverse shear stresses in composites [24]. This theory was 
used to develop a composite box beam model by McCarthy and Chattopadhyay [27-29]. 
In this model, each wall of the box beam is analyzed as a composite plate using a refined 
higher-order displacement field [30]. Continuity between the displacement fields is 
enforced at the four corners throughout the thickness of each plate. The model accurately 
captures the transverse shear stresses through the thickness of each wall while satisfying 
stress free boundary conditions on the inner and outer surfaces. The formulation 
approximates three-dimensional elasticity solution so that the beam cross-sectional 
properties are not reduced to one-dimensional parameters. Both inplane and out-of-plane 
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warping are automatically included in the formulation. The finite element method is used 
to solve the governing equations of motion. The model has been validated extensively 
for thin- and thick-walled composite laminates through comparisons with experimental 
results, other appropriate theories and three-dimensional finite element analysis using 
brick elements [27-29]. This method is used in the present work for accurate and 
efficient composite structural analysis of aircraft wings. 

1.1.3 Aeroelastic Analysis 

Aeroelastic analysis plays a vital role in the design of a high performance aircraft. 
History of the U.S. Supersonic Transport program shows that the entire aircraft design 
process was driven by the aeroelastic design cycle [31]. To effectively integrate 
aeroelastic analysis with the design of composite wings, computationally efficient yet 
analytically rigorous methods are necessary. The key issues associated with the 
aeroelastic stability analysis include structural dynamic calculations for natural frequency 
and mode shapes, unsteady aerodynamic computations to obtain generalized aerodynamic 
forces and flutter calculation methodology. These issues are briefly discussed next. 

The objective of the aeroelastic analysis is to identify the flight condition (velocity 
and atmospheric density or altitude) at which the aeroelastic system is neutrally stable. 
At this condition, the system is purely oscillatory and the aerodynamic loads calculated 
for simple harmonic motion are adequate. The technique for the prediction of three- 
dimensional unsteady aerodynamic forces for purely oscillatory motion is well 
developed. In the present research, the generalized aerodynamic forces are computed 
using the constant-pressure lifting surface method [32] at a given Mach number for 
specified values of reduced frequency k, assuming simple harmonic motion. This method 
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is based on linearized aerodynamic potential theory. The lifting surface is divided into 
small trapezoidal panels, with unknown constant pressure, arranged in strips parallel to 
the free stream. The downwash boundary condition calculated from and deflection of 
structural modes is satisfied at mid-span three-quarter chord point of each panel, leading 
to the determination of the unknown pressures and the aerodynamic forces. 

Assuming that the aeroelastic system performs a simple harmonic motion and the use 
of purely oscillatory air loads, leads to the classical V-g method of flutter prediction [33] 
which has been extensively used. However, the method is iterative in nature which 
reduces the efficiency and the results are accurate only at the flutter boundary. To gain an 
insight into the physical phenomena leading to flutter, it is necessary to obtain valid 
damping and frequency history. The Laplace domain method of flutter analysis [34] 
produces root-loci of the system which affords such an insight. The principal difficulty in 
implementing the Laplace domain method lies in obtaining the air loads for arbitrary 
motions, since aerodynamic calculations are well developed only for simple harmonic 
motions. This problem is overcome through the use of rational function approximations 
(RFA’s). The generalized aerodynamic forces contain transcendental terms when 
expressed as a function in the Laplace domain. To obtain a finite number of terms, the 
aerodynamic forces are approximated by a rational function of the nondimensionalized 
Laplace variable p. Several formulations of the RFA's are available in the literature [34- 
37]. The capabilities of these formulations have been extended and their performances 
compared in [38]. In the current work, the 'least-squares' approach of [35-36], which has 
been used by many researchers, has been adopted. 


1.1.4 Aerodynamic Analysis 
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To obtain aerodynamic efficiency (lift-to-drag ratio), it is necessary to compute lift 
and drag for the wing. Wing drag at high subsonic or transonic Mach numbers comprises 
induced drag, skin-friction drag and compressibility (wave) drag. The panel method 
based on the constant-pressure lifting surface [32] is also capable of predicting steady 
normal forces. The components of the normal force along the freestream and the 
direction perpendicular to it yield induced drag and lift, respectively. The skin-friction 
drag arises from the viscous effects in the mostly turbulent boundary layer adjacent to the 
wing surface. The turbulent boundary layer problem is a very difficult one to solve 
theoretically. Instead, empirical formulae developed by von Karman or Schlichting [39] 
are often employed to compute turbulent skin-friction drag. In the current study, the skin- 
friction drag is calculated using the Schlichting empirical formula, corrected to include 
the Mach number effect [40]. 

The compressibility drag refers to the pressure drag resulting from increase in Mach 
number above low subsonic value. At high subsonic or transonic Mach numbers, shocks 
develop on the top of the wing due to increased airflow velocity, which leads to 'drag 
rise'. The drag rise can be analytically predicted only through the use of sophisticated 
nonlinear computational aerodynamic analysis, since the linear analysis in this Mach 
number regime produces completely incorrect results [40]. Therefore, an empirical 
method described in [41] has been used in the current work to obtain the compressibility 
drag. This method, called the 'crest-critical Mach number method', has been used in other 
optimization studies [9, 42]. In this approach, the free stream Mach number which gives 
sonic flow at the highest point on the airfoil (tangent to the free stream) is first 
determined. This Mach number, called the crest-critical Mach number, is a function of 
airfoil thickness-to-chord ratio, lift coefficient and wing sweep. The compressibility drag 
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is then obtained using empirical relationship between the ratio of freestream Mach 
number to the crest-critical Mach number and drag rise. 

1.1.5 Smart Composite Wing Design 

Smart materials respond to an external stimulus, such as electric field, by changing 
their shape. When attached to a host structure, they cause deformation of the structure. 
The feasibility of using smart structures is increasing because of the availability of smart 
materials commercially, ease of integration with laminated structures, potential of large 
performance enhancement and advances in related fields [43]. Piezoelectric materials are 
popular for aeroelasti c/aerodynamic and vibration control [44]. When a piezoelectric 
material is stressed mechanically by a force, it generates an electric charge. Conversely, 
when an electric field is applied, the material elongates or shortens depending on the 
polarity of the applied electric field. A piezoelectric element is therefore capable of being 
used both as actuators and sensors. Currently, most piezoelectric devices utilize lead 
zirconate titanate (PZT), which is a piezoceramic material. The desirable properties of 
PZT include a high level of piezoelectric activity, a wide frequency range and first-order 
linearity between applied voltage and induced strain [43]. 

Composite wings with either embedded or surface bonded PZT actuators / sensors 
have been investigated by researchers. Heeg [45] demonstrated analytically and 

experimentally that piezoelectric materials can increase the flutter speed of a simple two 
degree of freedom wing model. Song et al. [46] showed that incorporation of 
piezoelectric layers in a wing can improve both divergence instability and aeroelastic lift 
distribution. Paige et al. [47] used piezoelectric actuation to control panel flutter. The 
structural modeling issues associated with piezoelectric actuation of composite plates 
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have also been considered. The displacement field based on the classical laminate theory 
is used by Chopra [43] and Crawley and Lazarus [48]. Chattopadhyay and Seeley [49] 
and Seeley [50] have used the higher-order theory (HSDT) for modeling composite plates 
with piezoelectric actuators/ sensors. The present research uses the above higher-order 
plate theory to model composite box beams with surface bonded piezoelectric actuators. 
The formulation of the box beam (with taper and sweep) is similar to that of McCarthy 
[27] and McCarthy and Chattopadhyay [28, 29]. The governing equations of motion are 
solved using the finite element method to be able to address realistic wing geometries. 

The placement and the number of actuators necessary for improved aeroelastic control 
require the use of formal optimization techniques. Several investigations have been 
reported which address the issue of actuator placement using both deterministic [51-52] 
and heuristic [53-55] approaches. For aeroelastic control, actuators placed at wing root 
have been shown to be most effective [52, 55]. Power consumption is also an important 
issue related to the piezoelectric actuation of structures, especially for active control with 
multiple actuators/sensors [56]. 

In the present research, a new procedure for the multidisciplinary design optimization 
of smart composite wings has been developed, incorporating optimization and analysis 
methods discussed above. The principal load carrying member of the wing is modeled as 
a composite box beam with surface bonded piezoelectric actuators. The optimization 
problem is formulated with the coupling of structural, aerodynamic, aeroelastic and 
control (passive) design criteria. The higher-order theory has been used for the wing 
structural analysis and its impact on aeroelastic results have been demonstrated through 
comparisons with those obtained using CLT. The effect of composite ply orientations on 
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flutter and divergence speeds has been studied The developed MDO procedures are 
demonstrated through applications to wing design. 

1.2 Multidisciplinary Design Optimization of Turbomachinery Blades 

High performance aircraft engine components operate under severe aerodynamic, 
thermal and structural environments. The design of the blade profile is one of the major 
aspects in engine design [57]. Engine performance is strongly affected by the 
aerodynamic efficiency of the blades, which can be enhanced through efficient design of 
the blade external shape. Sharp fluctuations in the blade (suction and pressure) surface 
Mach number can lead to flow separation resulting in loss of aerodynamic efficiency. 
Airflow velocity also impacts blade cooling and temperature distributions. The 
maximum and average temperatures of the blade are desired to be minimum, as the 
structural integrity and engine life are affected by these temperatures. From the structural 
point of view, it is important to maintain blade stresses and vibration levels within the 
allowable limits. Therefore, efficient blade design is a multidisciplinary problem that 
requires the integration of several disciplines such as aerodynamics, heat transfer and 
structures. 

The direct design method, in which the designer changes the blade geometry 
iteratively until desired performance is achieved, affords direct control of the blade design 
parameters. However, this method is very laborious and requires considerable insight 
[58]. In the inverse design method [59], the performance is specified in terms of velocity 
or pressure distributions to obtain the desired blade shape. This requires a knowledge of 
the desired velocity or pressure distribution. Also, the imposition of constraints is not 
easily applicable in inverse design. 
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Design optimization techniques have been used by many researchers to optimize 
airfoil shape for aircraft wings [60-61]. However, the application of optimization 
procedures to turbomachinery blade design has been rather limited. Chattopadhyay et al. 
[62] developed an optimization procedure for efficient design of turbine blades which 
successfully eliminated the leading edge velocity spikes while maintaining the tangential 
force coefficient. Aerodynamic analysis was performed using a two-dimensional panel 
code. The pressure and suction surfaces were approximated by polynomials, whereas 
circular and elliptic arcs were used to describe the leading edges. The procedure was 
further extended by Narayan et al. [63] to include heat transfer eriteria where coolant hole 
shapes and sizes were included as additional design variables. This multidisciplinary 
optimization procedure resulted in significant reduction in blade temperatures and smooth 
velocity distributions. For aerodynamic optimization, Goel et al. [64] used a combination 
of numerical optimization, hill climbing and genetic algorithm in an attempt to overcome 
the problem of local minima, since the turbine design problem is mutimodal. Turbine 
blade geometry was represented by Bezier-Bemstein polynomial [65]. Blade 
performance was measured by the distribution of the surface Mach number obtained 
through inviscid flow calculations. 

The current work presents the development and application of a new multidisciplinary 
optimization procedure incorporating more comprehensive analysis methods for the 
design of turbine blades. The procedure integrates aerodynamic and heat transfer design 
considerations, with mechanical constraints on blade geometry. Bezier-Bernstein 
polynomial is used to accurately represent airfoil shape with a relatively small number of 
design variables. The aerodynamic analysis is based on the thin shear layer 
approximation of the Navier-Stokes equations [66-67]. Grid generation is accomplished 
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by solving Poisson equations with arbitrarily specified inner and outer boundary points 
[68-69]. A finite element formulation is used to calculate blade interior temperatures. 
Total pressure and exit kinetic energy losses are minimized through a constrained 
multiobjective optimization formulation using the Kreisselmeier-Steinhauser (K-S) 
function approach discussed earlier. The maximum and average blade temperatures and 
certain geometric parameters of the blade are treated as constraints. 



2. Objectives 


The primary objective of the present research is to develop a multidisciplinary 
optimization procedure for the conceptual design of composite aircraft wings with surface 
bonded piezoelectric actuators. The optimization problem addressed involves the 
coupling of structural mechanics (including smart material), aeroelasticity and 
aerodynamics. The validity of the designs obtained using MDO procedures depends 
strongly upon the accuracy of the analytical methods used. Therefore, it is essential to 
integrate accurate and efficient analysis techniques to obtain meaningful optimum designs 
within a reasonable time. Since this multidisciplinary problem has multiple nonlinear 
objective functions and constraints, sophisticated optimization algorithm is required for 
solution. 

In the present research, the load carrying member of the wing is idealized and 
represented as a composite box beam. Each wall of the box beam is analyzed as a plate 
using a refined higher-order displacement field. This structural modeling accurately 
captures the transverse shear stresses through the thickness of each wall while satisfying 
stress free boundary conditions on the inner and outer surfaces of the beam. The present 
research extends the composite box beam model to include piezoelectric actuators bonded 
to top and bottom surfaces. Detailed structural modeling issues associated with 
piezoelectric actuation of composite structures are considered. The governing equations 
of motion are solved using the finite element method to analyze practical wing 
geometries. 

For the aeroelastic stability analysis, both the classical V-g method and the Laplace 
domain method are utilized. The V-g method gives accurate results at the flutter 
boundary, but requires iterative solution. The Laplace domain method involves 
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approximating generalized aerodynamic forces, but it produces root-loci of the system 
which give an insight into the physical phenomena leading to flutter/divergence. This 
method can be efficiently integrated within an optimization procedure. The steady and 
unsteady aerodynamic forces are obtained using linearized aerodynamic potential theory 
for compressible flows. 

The Kreisselmeier-Steinhauser (K-S) function technique is used to efficiently 
integrate the objective functions and constraints into a single envelop function for 
multiobjective optimization with continuous design variables. The resulting 
unconstrained problem is solved using the Broyden-Fletcher-Goldberg-Shanno algorithm 
for nonlinear optimization. The hybrid optimization method, though computationally 
expensive, includes continuous and discrete design variables simultaneously. 

The secondary objective of this work is to develop a multidisciplinary optimization 
procedure for the design of turbomachinery blades. Aerodynamic and heat transfer 
design objectives are integrated along with various mechanical constraints on the blade 
geometry. The blade geometry is represented by Bezier-Bemstein polynomials, which 
results in a relatively small number of design variables for the optimization. Thin shear 
layer approximation of the Navier-Stokes equation is used for the viscous blade-to-blade 
flow calculations. Grid generation is accomplished by solving Poisson equations. The 
maximum and average blade temperatures are obtained through a finite element analysis. 
Total pressure and exit kinetic energy losses are minimized, with constraints on blade 
temperatures and geometry. 


The specific goals of the current research are as follows; 



19 


1. Establish the significance of the refined higher-order displacement field on the 
aeroelastic stability of composite wings. Study the effect of composite ply 
orientations on flutter and divergence speeds. 

2. Extend the higher-order theory based composite box beam model to include 
piezoelectric actuators bonded to top and bottom surfaces, considering detailed 
structural modeling issues associated with piezoelectric actuation of composite 
structures. 

3. Develop a multidisciplinary optimization procedure for the conceptual design of 
composite aircraft wings incorporating accurate and efficient analysis methods 
and multiobjective optimization technique. The optimization problem is 
formulated with the objective of simultaneously minimizing wing weight and 
maximizing its aerodynamic efficiency. Design variables include composite ply 
orientations, ply thicknesses, wing sweep and piezoelectric actuator thickness. 
Constraints are placed on the flutter/divergence dynamic pressure, wing root 
stresses and the maximum electric field applied to the actuators. 

4. Develop of an accurate and computationally efficient optimization procedure for 
integrated aerodynamic and heat transfer design of turbomachinery blades. 



3. MDO Methodology for Smart Composite Wings 


The multidisciplinary design optimization of smart composite wings involves the 
coupling of structural mechanics (including smart material), aeroelasticity and 
aerodynamics. For the developed MDO procedure to be applicable to practical problems, 
the analysis and the optimization techniques must be computationally efficient and 
sufficiently rigorous. These methods are described in the following sections. 

3.1 Analysis 

3.1.1 Structural Modeling 

The load carrying member of the wing is represented as a single-celled rectangular 
box beam with taper and sweep (Figures 1 and 2). Piezoelectric actuators are bonded to 
top and bottom surfaces of the box beam. 

Y 



Figure 1. Wingplanform 
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The box beam is modeled using composite laminates to represent the four walls 
(Figure 3). The global coordinate system for the box beam is denoted (X, Y, Z) and the 
local coordinate system for the i-th wall is denoted (xi, yi, zi). The subscript 'i' is omitted 

for convenience in the rest of the dissertation. For each of the individual walls of the box 
beam, the inplane displacements are represented as cubic functions of the thickness 
coordinate and the transverse displacement is assumed constant through the laminate 
thickness. The higher-order displacement field [30] described in the local coordinate 
system (Figure 3) is as follows. 

u(x,y,z) = Uo(x,y) + zy^(x, y) + y) + 

v( X, y , z) = V o( X, y) + zvj; y (x, y ) + z^^y (x, y ) + z^^y (x, y) 

w(x,y) = Wo(x,y) (1) 

where Uq, Vq and Wq denote the displacements of a point (x, y) on the midplane and 
and v)/y represent the rotations of normals to the midplane about the y and x axes, 
respectively. The higher-order terms . Cx ’ represent beam warping in each 

plane. 
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3 PZT actuator 

(b) 


Figure 3. Smart box beam construction 

The displacement field must satisfy the conditions that the transverse shear stresses, 
and dy^, vanish on the plate top and bottom surfaces. 

<^xz(x,y,±-) = 0 
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ay^(x,y,±^) = 0 (2) 

where h is the plate thickness. For orthotropic composite plates, these conditions imply 
that the corresponding strains be zero on the surfaces. This yields the following relations. 


^x =0 


=0 

^ 4 SWn 


4 ,0Wn 

Cv ~ ~ M^v) 

3h^ ^ ^ 


( 3 ) 


The refined higher-order displacement field is now written as follows. 

4z^ ^awo 

u = Uo + z\]i^ 

4z^ , awp 

v=Vo + zv|;y-^(-^+M/y) 

W = Wq 


The functions Uq, Vq, Wq , v|^x v|/y represent unknown functions of x and y. It is 
important to note that the displacement field for the refined higher-order theory (Equation 
(4)) has the same number of dependent variables as those used in the first-order shear 
deformation theory (FSDT), although inplane displacements are cubic functions of the 
thickness coordinate. By making substitutions, 


Vx = ‘i'x 



dx 


Vy =<t>y - 


5wq 

dy 


( 5 ) 
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the following modified form of the displacement field is obtained. 



W = Wq 


where and (j)y represent the additional rotations due to shear deformation about the y 

and X axes, respectively. It should be noted that the displacement field for the classical 
laminate theory can be obtained by setting and (j)y , equal to zero in equation (6). 

Since there are only six unique values of the stresses and the strains due to symmetry, 
the following notation is used to define these quantities. 



where x, y and z correspond to 1, 2 and 3 directions, respectively. Assuming small 
displacements and rotations, a linear strain-displacement relationship is used. 
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^4 = 


du 

^ 2 =' 

5x 

dv 5w 
dz dy ’ 


dv 

dy'' 


— 


dw 

Si = — = 0 

dz 


5u 5w 
dz dx ’ 


du dv 
= + 
dy dx 


( 9 ) 


In the above equation, s, and £2 inplane normal strains, £3 is transverse normal 
strain, £4 and £5 are transverse shear strains and £g is inplane shear strain. Using the 
above strain-displacement relationships (Equation (9)) and the refined higher-order 
displacement field (Equation (4)), the strains can be expressed in terms of midplane 
displacements and curvatures as follows. 

£| = £]* + ZK® + 

0 , 0,32 

£2 = £2 + ZK2 + Z K2 

0,0 3 2 /I f\\ 
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The generalized Hooke's law is used to relate the stresses and the strains. For 
laminates made of orthotropic materials (with elastic symmetry parallel to the 1 -2 plane) 
the constitutive relation [70] is written as 

^1 Qii Qi 2 0 8i 

^ CT2 ’ = Qi2 Q22 0 ' 82 " 

0 Q66J1^^6- 


H Q44 0 1^4! 
_0 Qss-l^sJ 


(12) 


where Qy , the plane-stress reduced elastic constants in the material axes of the laminate, 
are related to the material engineering constants by the following equations. 

Qll= '^ 1 , Qj, = _S2— 

I-V12V21 l-''l2V2j (13) 

Q44=G23’ Q 55 “ G] 3 , Q66~G|2 


In the above equations, there are only five independent elastic constants. For laminates 
consisting of multiple plies at different orientations, it is convenient to use the 
transformed elastic coefficients [70] in the laminate coordinate system (x, y, z). After 
transformation to the laminate coordinate system, the constitutive relations can be written 
as 

O'] Qii Qi2 Qi 6 i| 

< q2 " = Qi2 Q22 Q26 ' ^2 ' 

'^6-1 _Q|6 Q26 Q66-^^6' 

1^4 1 _ Q44 Q45 1^4 1 

Kj Lq 45 QssJl'sJ (14) 


where Qy are the transformed elastic constants. 
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Figure 4. Composite laminate with surface bonded piezoelectric actuator 


The constitutive relation for any ply of piezoelectric material is written as follows 
[48], 

ai-Qy(ij-Aj) (i,j = l,2,6) (15) 

where A j is the induced strain due to piezoelectric actuation. Equation (15) is applicable 

if the piezoelectric material thickness is small compared to the plate thickness (Figure 4). 
The actuation strain vector Aj contains inplane normal and shear strain components, and 

can be treated similar to thermal strain in the elasticity formulation. The induced strain is 
used to control extension, bending or twisting of a laminate. Generally, an electric field 
is applied through the thickness of the PZT used as an actuator. Also, the PZT materials 
are isotropic and, therefore, its orientation has no effect on the material properties. 
Denoting the thickness direction as '3', the induced strain for PZT material is obtained 
from the following relation. 

^31 

' A2 ' = ' d 3 i 'E3 

Aei io ^ 


( 16 ) 








28 


In equation (16), E 3 is the applied electric field and d 3 , is the mechanical-electrical 
coupling coefficient of the material. The applied electric field is obtained by dividing the 
applied voltage by the thickness of the PZT layer. The maximum applied electric field 
must be smaller than the coercive field to avoid depoling the material. 

In the 'equivalent single layer' approach of composite analysis [70], the laminate 
stress resultants are obtained by integrating the stresses through the thickness as follows. 
The stresses also include effects due to piezoelectric actuation. 

h/^ . 

(Nj, Mj, Pj)= JaYl, z, z Uz (i = l,2,6) 

-h/2 

‘"'I / 9\ 

(Q2,R2)= {04(1, Z )iz 

-h/2 

(Q„R,)= laA x^z 

-h/2 ^ ^ (17) 

The first three terms (Nj, M j and Pj) are the inplane terms, which can be decomposed as 
follows. 

Nj =Nf -Nf 

Mi=Mf-Mf (i = l,2, 6) (18) 

p. = pA _ p.P 

The first terms on the right hand side equation (18) (superscript 'A') represent the stress 
resultants due to actual strain field, whereas the second terms (superscript 'P') represent 
the resultants from piezoelectric actuation. 

The stress resultants due to the actual strain field can be written in terms of the elastic 
constants, Qjj , and strains, Sj , as follows. 
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h/2_ 

(Nf , Mf , Pi^)= i Qiji|(l, z, z^)iz (i = 1 , 2 , 6) 


Substituting the expression for strains from equation ( 10 ), the plate constitutive relations 
are expressed as follows. 
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where the plate stiffnesses are defined as follows 

(Aij, By. Dy, Ey, Hy)= J Qy(l, z, z\ z^ )lz (i,j = 1,2,6) 

-h/2 

. vh/2 / 0 A\ 

(Ay,Dy,Fy)= ^ Qy (l, Z^ z'^)lz (i,j = 4 , 5 > 


The stress resultants due to piezoelectric actuation can be similarly defined as follows. 
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(Nf, Mf, Pn= TQijAi(U z, z^)iz (i = 1, 2, 6) (24) 

^ -h/2 

where the stiffness matrix, Qjj , includes the elastic constants of the piezoelectric material 
only. Using the stiffness matrices for the PZT material, the stress resultants can be 
expressed as 



where the PZT stiffness matrices are given by the following equation. 

h/2 

(Aij,Bjj. Ey)= I Qij(l, z, z^)lz (i,j= 1,2,6) (26) 

-h/2 

The box beam equations of motion are derived using the Hamilton’s principle [24] 
which assumes the following form, in the absence of any nonconservative forces. 

6l[K-(V+U)]dt = 0 (27) 

t| 

where 5 represents the variation and K is the kinetic energy, V is the potential energy due 
to external forces, U is the strain energy and t] and t 2 are initial and final times. 
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respectively. Using variational principles, equation (27) may be written in terms of the 
individual wall (plate) quantities as follows. 

Uf N 1 


J 


I5Kj -5Vj -5Ui dt = 0 


(28) 


l|Li=l 


where N is the total number of walls (N=4 for a box beam). The variations of the elastic 
strain energy, the potential energy and the kinetic energy, in each plate, are written as 
follows. 

h/2 

5U = j {(aibEj + a26s2 +<74684 +a56e5 +a58e6)iAdz 
-h/2A 
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sh/2 

6K = — J I pru + w + MAdz (29) 

2 -h/2A ^ ^ 


where (■) denotes differentiation with respect to time, A is the plate area, p is the 
material density and p(x,y) is the distributed load. Substitution of equation (29) into 
equation (28), integration and collection of the coefficients of various terms, yields the 
following equations of motion. 
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,(30) 


where the stress resultants Nj, Mj, Pj, Q; and R^ are as defined in equation (18) and the 
inertia terms are defined as 



For the solution of the equations of motion, appropriate natural or essential boundary 
conditions [70] needs to be specified. 

The finite element method is well suited to solve the above equations of motion, 
accounting for discontinuities in the material properties and complex geometry. The 
discretized equations of motion are obtained using a two-dimensional finite element 
formulation in the local coordinate system of each individual wall. To maintain 
continuity of displacements and to ensure that the walls remain normal to each other after 
deformation, appropriate constraints are imposed on displacements and rotations of 
individual walls at the four comers of the box beam cross section (Figure 3). A four 
noded plate element is used with 1 1 degrees of freedom per node. This element is C ^ 
continuous in the zeroth order displacements (uq, Vq and Wq ) and continuous in the 
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higher-order terms (\|/x and v|/y). The nodal degree of freedom vector is defined as 
follows. 


q = 


5uq 5U| 


•* 0 ’ 


5x ’ 5y ’ 5x ’ 5y ’ dx ' dy 


5vq dv 
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Wo 


5wn 5w, 


IT 




(32) 


Using shape functions, each of the unknowns is interpolated over the element as follows. 
U®(x,y) = ZN-(x,y)qf (33) 

i=l 

where n is the number of nodes (n=4), Nf is the shape function and superscript 'e' denotes 
elemental representation. Bilinear shape functions are used for Uq, Vq, v|;^ and \\iy, 
whereas a 1 2 term cubic polynomial is used for the transverse displacement Wq . The 
variations of kinetic, strain and potential energies, in terms of the elemental quantities, are 
obtained and the elemental matrices are assembled leading to the following equations of 
motion. 


m|+ Kq = f + Fp 


where the global mass and stiffness matrices are given by 
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M = E 

i=l 


jpS'SdV 

LV 


_ 4 
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i=l 


jB'QBdV 


(34) 


(35) 

(36) 


and the force vectors are as follows. 


f = zl 

i=t 


jNVx,y)dA 
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( 37 ) 
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4 

Fp =Z 


Jn’^n’^cIA 


La 


(38) 


In equations (35) to (38), the summation is performed over walls 1 through 4 (number of 
walls=4 for a box beam) and V and A represent volume and surface areas; respectively. 
The material density is denoted p and p(x,y) is the distributed load (such as, aerodynamic 

pressure). The matrix Q is the material stiffness matrix and matrices B and S relate the 

T 

nodal degree of freedom vector to strains and displacements, respectively. The term N 
denotes transpose of the shape functions, N*’ is the stress resultant due to piezoelectric 
actuators and Fp is the corresponding force vector. Since the accurate prediction of 
structural damping is difficult, it is ignored in the present analysis. Also, its effect is 
usually small on aeroelastic stability and ignoring it generally leads to conservative 
results [34]. 

3.1.2 Aeroelastic Analysis 

The objective of the aeroelastic analysis is to obtain velocity and atmospheric density 
or altitude where the aeroelastic system is neutrally stable. The V-g (velocity - damping) 
method of flutter prediction [33] is the classical method which has been extensively used 
over the past decades. It assumes that the aeroelastic system performs a simple harmonic 
motion and uses the purely oscillatory air loads. Since these assumptions are valid at the 
flutter condition, this method yields accurate results for flutter boundary. 


V-g method 


The equations of motion for the wing structural-aerodynamic system is solved with 
the aeroelastic forces represented by the force vector f . Assuming simple harmonic 
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motion, wing elastic deformation and aeroelastic forces are represented by q = qe‘“' and 
f = fe'“* , respectively. Substituting the above into equation (34) yields 

(-(O^M+K)q = f (39) 

where the stiffness term, K, now also includes the piezoelectric actuation effects. The 
unsteady aerodynamic force, f , can be expressed as a linear combination of q as follows. 

f = q^F(ico)q (40) 

where F(ico) is the aerodynamic influence coefficient and denotes freestream 
dynamic pressure. Substituting for f in equation (40) results in the following equation. 

(-o)^M+K-q^F(ico))q = 0 (41) 


Equation (41) represents an eigenvalue problem and the solution of the following 
determinant 


-co^M+ K-qj^F(ico) 


= 0 


(42) 


provides the roots which determine the stability of the system. To solve the above 
problem, artificial damping, g , is introduced and equation (42) is rewritten as 

-co^M + (1 + ig)K - q^F(i 03 )|= 0 (43) 


The problem size for flutter solution can be reduced by using the modal approach. It 
is well known that only the low frequency modes govern wing flutter characteristics and 
the high frequency modes have little effect on the flutter solution [33]. Using the first 
several low frequency modes, equation (43) is transformed into modal coordinates and is 
rewritten as follows. 
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|-(oW+(l + ig)K*-q^F*(ico)|=0 (44) 

where M , K and F are the generalized mass, stiffness and aerodynamic influence 
coefficient, respectively. The generalized aerodynamic forces are computed using the 
constant-pressure lifting surface method [32] at a given Mach number for specified values 
of the reduced frequency k, assuming simple harmonic motion. This method is based on 
the linearized aerodynamic potential theory [33]. The lifting surface is divided into small 
trapezoidal panels, with unknown constant pressure, arranged in strips parallel to the 
freestream. The downwash boundary condition calculated from the deflection of 
structural modes is satisfied at mid-span three-quarter chord point of each panel. The 
procedure is implemented in the 'ZONA6' computer code [71] which is used in the 
present research. The solution of equation (44) yields the variations of damping and 
frequency with respect to the freestream velocity. At the flutter point, the artificial 
damping, g=0. 

Laplace domain method 

The V-g method gives accurate results at the flutter point, but it does not generate 
reliable frequency and damping histories since the assumed harmonic motion is not valid 
at other conditions. Also, this method requires iterative calculations to arrive at the 
'matched flutter point', which poses difficulty in using it within an optimization 
procedure. The Laplace domain method of flutter analysis [34] produces valid damping 
and frequency history and, thereby, affords an insight into the physical phenomena 
leading to flutter. This method is non-iterative and suitable for automated optimization 
procedures. However, the Laplace domain method requires the air loads for arbitrary 
motions. Since aerodynamic calculations are well developed only for simple harmonic 
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motions, approximations are necessary to define air loads. These approximations reduce 
the accuracy of flutter prediction, but the method is still very useful. 

Laplace transformation of the governing equation of motion (Equation (34)) yields 
the following. 

(s^M + K)q(s) = f(s) (45) 

where s denotes the complex Laplace variable. The aeroelastic load f(s) can be expressed 
as a linear combination of q(s) as follows. 

f(s) = q^F(p)q(s) (46) 

where p is the nondimensionalized Laplace variable and F(p) can be regarded as the 
aerodynamic transfer function. Substituting for q(s) in equation (45) gives 

(s"M + K-q^F(p))q(s) = 0 (47) 

Equation (47) can be transformed into modal coordinates and rewritten, using several low 
frequency modes as follows. 

(s2M%K*-q^F*(p))^(s)=0 (48) 

where M*,K* and F* are the generalized mass, stiffness and aerodynamic forces 
respectively and ^(s) denotes the generalized coordinate. 

The unsteady aerodynamic forces are obtained assuming simple harmonic motion [32, 
71] similar to the V-g flutter calculations. However, F*(p) contains transcendental terms 
when expressed as a function in the Laplace domain. To obtain a finite number of terms, 
the aerodynamic forces are approximated by a rational function of the 
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nondimensionalized Laplace variable p. Following the method of [35-36], the 
approximating function for F (p) is expressed as 

F(p) = Aq + A]P + A2P^ + Z(A|^_2) (49) 

1=1 P+Di 

where n is the number of partial fractions (order of fit). The partial fractions approximate 
the time delays inherent in unsteady aerodynamics [36]. The denominator coefficients b| 
in equation (49) are selected from the range of reduced frequencies for which unsteady 
aerodynamic forces are computed. Substituting ik for p along the imaginary axis, the 
coefficients, Aq, A|, etc. are computed from the available generalized forces sueh that 
the approximation error is minimized in the least-squares sense. Substituting equation 
(49) into equation (48), the equations of motion are reduced to a series of 6N first-order 
equations [36] of the form 

sX = AX (50) 

where N is the number of normal modes used. The eigenvalues of matrix A provide the 
roots of the flutter equation. For stability of the system, all real roots should be negative. 
At flutter condition, one of the roots is purely imaginary. The above procedure is 
implemented in the computer code 'Interaction Structures, Aerodynamics, and Controls' 
(I SAC) [72] and is used for flutter analysis in this research. 

3.1.3 Aerodynamic Analysis 

Lift and drag for the wing are calculated to define its aerodynamic efficiency (lift-to- 
drag ratio). The constant-pressure lifting surface method [32, 71], used for unsteady 
aerodynamic calculations, can also predict steady normal forces with the value of the 
reduced frequency being set to zero. Induced drag and lift are obtained by resolving the 



39 


normal force along the freestream and the direction perpendicular to it. The other two 
components of wing drag at high subsonic or transonic Mach numbers include skin- 
friction drag and compressibility (wave) drag. 


The skin-friction drag arises from the viscous effects in the mostly turbulent boundary 
layer adjacent to the wing surface. First, a flat-plate skin-friction drag coefficient (Cf) is 
calculated, which is multiplied by a 'form factor' (FF) to account for the viscous 
separation effects. The flat-plate skin-friction coefficient depends upon the Reynolds 
number, the Mach number and the skin roughness. For turbulent flow, which is generally 
the case at high subsonic or transonic Mach numbers, the flat-plate skin-friction 
coefficient is determined using the Schlichting empirical formula, corrected to include the 
Mach number effect [40], as follows. 


0.455 

(logioR)^^*(l+0.144M2)®-^^ 


(51) 


where M is the Mach number and the Reynolds number (R) is defined by 

R = pVf/|a (52) 

where p is the atmospheric density, V is the velocity, f is the mean aerodynamic chord 
length and p is the coefficient of viscosity. For relatively rough surfaces, the friction 
coefficient is higher. Therefore, a 'cut-off Reynolds number' is determined using the 
mean aerodynamic chord, f , and a surface roughness parameter, k . 

Rcut-off = 44.62(f/K)'°^^M'-'^ (53) 

The lower of the cut-off Reynolds number and the actual Reynolds number is used in 
equation (51). The form factor for the wing is obtained as follows. 
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FF = 



(54) 


where (x/c)^^ denotes the chordwise location of the airfoil maximum thickness point, 
f-1 is the airfoil thickness to chord ratio and refers to sweep of the maximum 


vc/ 

thickness line. The skin-friction drag coefficient can now be calculated from the 
following equation. 
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(55) 


where 


^vvet 

Sref 


1.977 + 0.52 




>0.05 


(56) 


represents the ratio of 'wetted area' to 'reference area'. 

The compressibility drag refers to the pressure drag resulting from increase in Mach 
number above low subsonic value. At high subsonic (or higher) Mach numbers, the local 
velocities on the upper surface of the wing may become sonic or even supersonic. This 
may lead to shock formation on the top of the wing which increases drag due to reduction 
in the total pressure through shock waves. Drag may also increase due to thickening or 
separation of the boundary layer as a result of the severe adverse pressure gradient caused 
by the presence of shocks. An empirical method, called the 'crest-critical Mach number 
method' [41], has been used in this research to obtain the compressibility drag. The 'crest' 
is the point on the airfoil upper surface to which the freestream is tangent. In this 
method, the freestream Mach number which gives sonic flow at or behind the crest of the 
airfoil is first determined. The crest-critical Mach number is a function of airfoil 
thickness-to-chord ratio, lift coefficient and wing sweep. At a Mach number 2-4 percent 
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(depending on wing sweepback angle) higher than the crest-critical Mach number, the 
drag rises abruptly. This Mach number is called the 'drag divergence Mach number', 
which can be determined from the crest-critical Mach number. The compressibility drag 
coefficient is then calculated using the ratio of freestream Mach number to the crest- 
critical Mach number and empirical data (based on existing transport aircraft) for the 
increment in drag. 

3.2 Optimization Technique 

The optimization techniques described next are suitable for multiobjective 
optimization with objective function/constraint approximation and continuous/discrete 
design variables. 

3.2.1 K-S Function Approach 

The K-S function technique [16] is used to efficiently integrate the objective function 
and constraints into a single envelop function. The resulting unconstrained nonlinear 
optimization problem is solved using the Broyden-Fletcher-Goldberg-Shanno (BFGS) 
algorithm [14]. The derivatives of the objective functions and the constraints, with 
respect to the design variables, are calculated using the forward finite difference 
technique. 

Using the K-S approach, each of the original objective functions is transformed into 
reduced objective functions as follows. 

» F (0>) 

F;(0) = -^-1.0-g^,,<0, 


i=l,...,NOBJ 


(57) 
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where F® is the value of the original objective function Fj calculated at the beginning of 
the cycle, O is the design variable vector and g^^^ is the value of the largest constraint 
and is held constant during each cycle. It is assumed that the constraints gj(0) (j = 1,2, 

, NC) are scaled to lie between -1 and 1. The reduced objective functions are 
analogous to constraints. A new constraint vector fm(d’) (m=l,.., M where M = NC + 
NOBJ) is defined by combining the original constraints and the reduced objective 
functions. The new objective function is defined as follows. 

1 M 

Fks(®) = fm.x +-loge Ze”"'- ' (58) 

P m=l 

where f„,ax is the maximum value of fm(d>) and remains constant during the 
optimization cycle. The scalar multiplier p is similar to draw-down factor of penalty 
function methods. Larger values of p move the K-S function curve closer to the largest 
constraint [14]. The initial value of this parameter is user supplied and its value is 
increased as the optimization proceeds. The new objective function Fj;g(0), which 
represents an envelope function representing the original objective functions and 
constraints, can now be minimized using any unconstrained optimization technique. In 
this research, the search direction vector is obtained using the Broyden-Fletcher- 
Goldberg-Shanno (BFGS) algorithm which is a quasi-Newton technique [14]. This is 
followed by a one-dimensional search for the minimum of the K-S function using the 
three point quadratic approximation for step length calculation [14]. 

3.2.2 Approximation Technique 

During the one-dimensional search to minimize the composite K-S function, several 
evaluations of the objective functions and constraints are necessary. It is computationally 
expensive to carry out exact analysis all the time. Therefore, an approximation technique 



43 


is used to provide objective function and constraint values during the one-dimensional 
minimization. The two-point approximation technique [18], which has been found to be 
well suited for nonlinear optimization problems, is applied. This technique uses the 
gradient information from the previous and current design cycles to obtain the exponent 
used in the expansion. The technique is formulated as follows. 


NDV 

F(0)=F(Oo)+ S 

n=l 


f A V" 


1.0 




Pn 


(«>o) 


(59) 


where F(0) is the approximation to the objective function F(0) at a neighboring design 
point <1) (vector of <])„ design variables) based on its values and gradients at the current 
design point Oq and the previous design point O, . The exponent p„, considered as a 
'goodness of fif parameter, explicitly determines the trade-offs between traditional and 
reciprocal Taylor series expansions. In the limiting case of p„ = 1, the expansion is the 
first order Taylor series, and when p„ = -1, the two-point exponential approximation 
reduces to the reciprocal expansion form. The exponent p„ is obtained from the 
following equation. 


loge] 

la|>„ J 

’-logcj 

^ 5F 

(<t«o)} 

logel 
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>• 
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(60) 


The exponent p„ is defined to lie between +1 and -1. If any singularity is encountered, 
the exponent is set to +1 to obtain linear Taylor series expansion. A similar 
approximation is obtained for the constraint vector too. 


3.2.3 Hybrid Optimization 
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The inclusion of both continuous and discrete design variables significantly 
complicates the optimization problem. This is because the discrete design variables are 
not compatible with traditional gradient based optimization methods. Similarly, the 
continuous variables are not compatible with combinatorial optimization methods, such 
as branch and bound techniques, which require discrete values to operate. Therefore, a 
hybrid optimization technique developed by Seeley and Chattopadhyay [19] which 
combines both types of design variables is used. 

The general continuous/discrete optimization problem can be mathematically stated 
as follows. 

Minimize f(0(.,0<j) (61) 

Subject to g(Oj.,Oj). <0 (j = 1, 2,..., NCON) (62) 

Side constraints ^ 

(63) 

where f is the objective function, gj are the constraints, are the continuous design 
variables and Oj are the discrete design variables which can be selected from among a 
set of q preselected values. The hybrid optimization procedure is based on Simulated 
Annealing (SA) [73] where the design space is sampled by repeatedly perturbing the 
discrete design variables. At each iteration of the SA procedure, the objective function is 
minimized with respect to the continuous variables using a BFGS search algorithm [14]. 
This significantly improves the efficiency of the hybrid algorithm by directing the search 
using the gradient information when available. The constrained problem is formulated 
using a penalty function approach [14]. 



4. MDO Applications for Smart Composite Wings 


The analysis techniques and optimization methods described in the preceding chapter 
have been applied to the design of a composite wing, with and without smart materials 
[74-80]. The aeroelastic stability of composite wings was investigated using the V-g 
method in [74-75] and the Laplace domain method in [76], A hybrid optimization 
technique (combining both continuous and discrete design variables) was adopted for the 
wing designs investigated in [77-78]. Multidisciplinary optimizations of business jet type 
composite wing and smart composite wing are reported in [79] and [80], respectively. 
These applications, from the current research, are presented in the following sections. 

4.1 Aeroelastic Stability Using Higher-Order Laminate Theory 

First, it is necessary to investigate the effect of higher-order theory on aeroelastic 
stability of composite wings. To demonstrate this, aeroelastic analysis is performed for a 
simple scaled wing model shown in Figure 5. The wing semi-span and root chord are 90 
and 20 inch, respectively. The aspect ratio and the taper ratio of the wing are 12 and 0.5, 
respectively, and the mid-chord is unswept. It is assumed that the box beam is the 
principal load-carrying member of the wing which extends through the entire semi-span 
and is fixed at the root. The width and the height of the box beam are assumed to be 50 
and 10 percent of the local wing chord; respectively. The center of the beam coincides 
with the mid-chord of the wing. Each wall of the box beam is made up of eight 
unidirectional composite laminates (each consisting of several plies with identical 
orientation), which are symmetric about the mid-surface. The composite material (Gr-Ep 
T300) properties are listed in Table 1 [81]. The material density is multiplied by a factor 
of eight to account for the non-structural mass of the wing. The walls have a uniform 
spanwise thickness of 0.80 inch. The top and the bottom walls have (0/90/30/30)s lay-up, 
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while the vertical wall ply angles are (45/-45/45/-45)s- The direction of positive ply 
angle is indicated in Figure 5. 

Y 



Table 1. Material properties (Gr-Ep T300) 


Efi (msi) 

19.00 

E7(msi) 

1.50 

GLT(msi) 

1.00 

G'p'p(msi) 

0.90 

''lt 

0.22 

p (slugs/in^) 

1.68*10-3 


The structural analysis of the box beam is performed with 15 elements spanwise and a 
single element chordwise. The structural degrees of freedom for the box beam are 664 
and 844 using the CLT and the present approach, respectively. The first ten normal 
modes of vibration are used in the subsequent analysis and the natural frequencies are 
presented in Table 2. Although, bending-torsion coupling exists due to the ply 
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orientations used, for identification purposes only the dominant mode shapes are 
indicated in Table 2. Table 2 shows that the natural frequencies of the beam bending 
modes (Bl, B2, B3, B4 and B5) are somewhat higher for the CLT compared to the 
present approach. The torsional mode (Tl) exhibits the largest difference in the natural 
frequency (about three hertz), the CLT value being higher again. This is due to the fact 
that the transverse shear stresses are ignored by the CLT, but are efficiently modeled in 
the present approach. The frequencies of the chordwise bending modes (Cl, C2, C3 and 
C4), which are not used in the flutter analysis, are identical using both theories. 

Table 2. Natural frequency (Hz) and modes 


Mode 

Number 

Mode 

Type 

Higher-order 

Theory 

Classical 

Theory 

1 

Bl 

5.34 

5.37 

2 

Cl 

18.17 

18.17 

3 

B2 

27.75 

27.94 

4 

B3 

71.42 

71.99 

5 

C2 

84.07 

84.07 

6 

B4 

130.71 

132.36 

7 

Tl 

134.43 

137.23 

8 

C3 

202.32 

202.32 

9 

B5 

205.75 

207.70 

10 

C4 

256.42 

256.42 


Legend: B- Beam bending, C- Chordwise bending, T- Torsion 


Mode shapes for the first beam bending (Mode 1 ) and the first torsion (Mode 7) are 
presented in Figures 6 and 7 for the CLT and Figures 8 and 9 for the present theory. The 
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mode shapes are normalized such that the generalized mass equals unity. Mode 1 is a 
pure beam bending mode, whereas Mode 7 is a coupled fourth bending/first torsion mode 
in both cases. These two modes are selected for comparison because of their contribution 
in wing flutter as will be shown later. The response of Mode 1 is nearly identical (as are 
the natural frequencies, see Table 2) using both theories, but differences are observed in 
Mode 7. The present approach shows somewhat higher bending displacements in this 
mode, whereas the torsional displacements are slightly lower compared to CLT. 



Spanwise distance, in 


Figure 6. Bending mode shape (CLT) 



Spanwise distance, in 



Torsional mode shape Bending mode shape 


Figure 7. Torsional mode shape (CLT) 



Figure 8. Bending mode shape (present theory) 



Figure 9. Torsional mode shape (present theory) 
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The generalized aerodynamic forces are computed at M=0.95 with the wing divided 
into 48 panels. Flutter solution is obtained using the V-g method for various values of 
the atmospheric density (an input parameter) till the flutter speed and the air speed match 
each other ('matched flutter point'). The flutter results, in terms of variations of frequency 
and damping with airspeed, are presented in Figures 10 and 11 using both CLT and the 
present approach. 
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Figure 10. Frequency vs. Airspeed 

Flutter is often caused by the coupled bending and torsional motion of aircraft wings, 
wherein the frequencies of these two modes come close or coalesce around the critical 
flutter speed and the damping of either of the modes goes to zero [82]. Figure 10 shows 
the tendency of the first bending (Mode 1) and the first torsion (Mode 7) modes to get 
close to each other. For the present approach, the separation between the two frequencies 
is smaller and the tendency to come close is more pronounced compared to the CLT. 
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Figure 1 1 . Damping vs. Airspeed 


As a result, the damping of one of the modes (Mode 1) becomes zero at lower airspeed, as 
predicted by the present approach (Figure 11). The flutter speeds, using the present 
approach and the CLT, are 455 and 503 KEAS, respectively. This significant difference 
in flutter speed of 48 KEAS (about 10 percent) is very critical. Further, it is important to 
note that the CLT results are nonconservative. This establishes the need for using refined 
structural modeling techniques since flutter is a catastrophic phenomenon. 

4.2 Effects of Ply Orientations on Aeroelastic Stability 

To study the effects of elastic couplings on aeroelastic stabilty of composite wings, 
various ply orientations are investigated. The wing geometry is the same as in the 
previous section (Figure 5), with the exception that the walls now have an uniform 
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spanwise thickness of 0.16 inch and non-structural masses are not included in the 
analysis. 

Three types of ply lay-ups, (0/90/0/90)s, (0/90/30/30)s and (0/90/-30/-30)s, are used 

for the top and the bottom flanges of the wing. The vertical webs have the cross-ply 
orientations of (0/90/0/90)s in all cases. For the cross-ply lay-up there is no elastic 
coupling, whereas the (0/90/30/30)s orientation displays negative bending-twist coupling, 
that is, upward bending causes nose down twist. The stacking sequence (0/90/-30/-30)s 
displays positive bending-twist coupling, that is, upward bending generates nose up twist. 
The structural analysis of the box beam is performed with 10 elements spanwise and a 
single element chordwise. The first six normal modes of vibration are used in the 
subsequent analysis. The natural frequencies and the associated modes are listed in Table 
3. The first, third and fourth modes are beam bending modes and the fifth mode 
represents the first torsion mode in all cases. The second and the sixth modes are 
chordwise bending modes for (0/90/30/30)s and (0/90/-30/-30)s orientations. For the 

cross-ply case, the second mode is a chordwise bending mode, whereas the sixth mode is 
a torsion mode. 

Table 3. Natural frequency (rad/sec) and modes for different ply lay-ups 


Mode 

(0/90/0/90)s 

(0/90/-30/-30)s 

(0/90/30/30)s 

1 

75, B 

66, B 

66, B 

2 

289, C 

269, C 

269, C 

3 

331, B 

309, B 

309, B 

4 

761, B 

738, B 

738, B 

5 

875, T 

1028, T 

1028, T 

6 

1134,T 

1195,C 

1195,C 
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Legend: B- Beam bending, T- Torsion, C- Chordwise bending 

The generalized aerodynamic forces are computed at M=0.8 for values of the reduced 
frequency between 0 and 1 . The wing is divided into 48 panels. The Laplace domain 
method is used in the flutter calculations. Rational function approximation of the 
generalized aerodynamic forces (tabular data) is performed using four denominator 
coefficients. These coefficients are selected to be 0.2, 0.4, 0.6 and 0.8 and produce a very 
good fit. Figures 12-15 show a comparison of the tabular data with s-plane fit for real 
and imaginary parts of the two of the aerodynamic influence coefficients, F(l,l) and 
F(2,2). These comparisons are typical of the rest of the coefficients. The approximation 
error for the aerodynamic forces is generally less than one percent and it never exceeds 
six percent. Total least-square error for the rational function approximation equals 0.36 
percent. 



Figure 12. RFA of real F(l,l) aerodynamic coefficient 
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Figure 13. RFA of imaginary F(l,l) aerodynamic coefficient 



Figure 14. RPA of real F(2,2) aerodynamic coefficient 




Figure 15. RFA of imaginary F(2,2) aerodynamic coefficient 

The eigenvalues of the characteristic equation vary with dynamic pressure. A root 
locus constructed by varying the altitude for a given Mach number gives flutter dynamic 
pressure when one of the real roots becomes zero and the imaginary root is non-zero. 
Divergence is indicated when both the real and the imaginary roots of any one mode 
reduce to zero. The variations of frequency and damping with dynamic pressure provide 
an insight into flutter and divergence onset. Figures 16-24 show the dynamic pressure 
root locus, frequency history and damping history for the three different ply lay-ups. 
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Figure 16. Dynamic pressure root locus for (0/90/0/90)s ply angle 
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Figure 17. Frequency versus dynamic pressure for (0/90/0/90)s ply angle 
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Dynamic pressure, psi 

Figure 18. Damping versus dynamic pressure for (0/90/0/90)s ply angle 
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Figure 19. Dynamic pressure root locus for (0/90/-30/-30)s ply angle 
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Figure 20. Frequency versus dynamic pressure for (0/90/-30/-30)s ply angle 
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Figure 21 . Damping versus dynamic pressure for (0/90/-30/-30)s ply angle 
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Figure 22. Dynamic pressure root locus for (0/90/30/30)s ply angle 
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Figure 23. Frequency versus dynamic pressure for (0/90/30/30)s ply angle 
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Figure 24. Damping versus dynamic pressure for (0/90/30/30)s ply angle 

The cross-ply lay-up wing with (0/90/0/90)s ply orientation flutters at 4.88 psi when 
the frequencies of the first and the fifth modes come close and the damping of the first 
mode becomes zero (Figures 16-18). This lay-up has uncoupled elastic modes and the 
fifth mode is the first torsion mode for this wing. The wing with (0/90/-30/-30)s ply 

orientation does not show flutter, but due to a positive bending-twist coupling, it is prone 
to divergence. As the flight dynamic pressure increases, the frequency and the damping 
of the first bending mode reduce to zero at 2.66 psi (Figures 19-21). For the wing with 
(0/90/30/30)s ply angles, flutter or divergence is not encountered at the Mach number 

investigated (Figures 22-24). As mentioned earlier, this lay-up has negative bending- 
twist coupling which prevents divergence. Similar trends for the effects of ply orientation 
on flutter and divergence have been observed by other researchers [83]. 
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4.3 Application of Hybrid Optimization Technique 

The hybrid optimization technique has been applied in the continuous/discrete design 
variable optimization of a composite wing. The wing weight is used as objective function 
which is minimized with constraints on flutter/divergence speed and stresses at the wing 
root due to the specified air loads. Wing root chord and wall thickness are used as 
continuous design variables, whereas the ply orientations are treated as discrete variables. 


4.3.1 Optimization Problem 


The wing geometry, for this optimization problem, has unswept mid-chord (similar to 
Figure 5), an aspect ratio of 20 and a taper ratio of 0.5. The root chord is varied during 
the optimization and consequently, the wing span and area also vary. The objective is to 
minimize the weight of the box beam which represents the structural member in the wing. 
Therefore, in the optimization problem. 


f= W 


(64) 


where W is the weight of the box beam. The weight of the remaining components in the 
wing, such as the skin, is not considered in this work. Constraints are placed on the 
flutter speed and the maximum allowable stresses. The flutter/divergence speed (Vf) is 

constrained to be greater than 450 knots equivalent air speed (KEAS) at a flight condition 
of Mach=0.7 at sea level. This constraint is expressed as follows. 



(65) 
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The Tsai-Wu failure criterion [84] is imposed on the critical ply stresses at the root 
section where material failure is most likely to occur. The criterion states that failure 
initiates when the following inequality is violated. 

Fj. = FiGi +Fyajaj < 1 (i,j = 1,2,...,6) (66) 

where Oj represents stresses in the material coordinate system and Fi and Fij are related 
to the tensile and compressive yield strengths of the material (Table 4). This constraint is 
expressed as follows. 

g2=F,-l<0 (67) 

The static lift and drag are computed empirically [40] for the above flight condition at an 
angle of attack of 3.5°. This results in a lift force of 2350 lb and a drag force of 1 10 lb for 
the reference wing. This load is assumed to remain constant during the optimization. 
The spanwise load distribution is assumed to be elliptical over the wing planform. The 
root chord (cr) is defined as a continuous design variable as follows. 

c, =(j), 15"<(|), <25" (68) 

where the upper and lower bounds are also indicated above. The actual dimensions of the 
wing, such as the semi-span (s), the tip chord (ct) and the wing area (A) are all computed 
from Cr so as to retain the aspect ratio and the taper ratio values mentioned earlier. The 
box beam is constructed from Carbon-PEEK composite material [81]. Material 
properties are presented in Table 4. 
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The box beam walls are assumed to be made of several unidirectional composite plies 
with thickness tp=0.005”. The quantity N, which denotes the number of two-ply layers in 

each laminate, is an integer value and is defined as a discrete variable 

N = (j)2 (69) 

where (|)2 = [4,5,6,.., 14], The box beam wall thickness, h, is then determined as 

h = 2Ntp (70) 


where it is assumed that all four walls have the same thickness. Therefore, the total 
number of plies in each wall is 8 < 2N < 28 and h ranges between values of 0.04” < h < 
0.14”. 


Table 4. Material properties (Carbon-PEEK) 


El (msi) 

19.40 

Ex(msi) 

1.29 

GLT(msi) 

0.74 

GTT(msi) 

0.50 

gLT 

0.28 

density 

(slug/inch^) 

1.8x10-3 

Ultimate Strengths 


XT(ksi) 

309.0 

XC(ksi) 

160.0 

YT(ksi) 

11.6 

YC(ksi) 

29.0 

S(ksi) 

23.2 
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Four design variables (<l) 3 ,<t) 4 ,<l) 5 ,<t) 6 ) are used to describe the stacking sequence in 
each of the box beam walls as follows. 

Topwall [Qi/QiJn =[<f*3/<l’4]N 

Bottomwall [9i /02 ]n ~ 

Leftside wall [0i /02 ]n = [<t> 5 /<t> 6 ]N 

Right side wall [9i /02 ]n = <I>6]n 

where the material axes for the plies in each wall are referenced to their respective local 
coordinate systems (Figure 3). In the global coordinate system, the top and the bottom 
walls have the same ply angles, as is the case for the two side walls. The ply angles are 
allowed to vary from within a range of pre-selected values (|)j=[-45°, 0°, 45°, 90°]. All 
four walls are assumed to have the same thickness, which does not vary along the span or 
the chord. 

4.3.2 Results and Discussions 

Results obtained using the hybrid optimization procedure are presented in Table 5 and 
in Figures 25-28. Optimization results are compared with a reference design, which is 
selected based on engineering judgment. It should be noted that the optimum design is 
independent of the initial design due to the probabilistic nature of the hybrid optimization 
procedure. The penalty function optimization iteration history is presented in Figure 25 
at each iteration of the simulated annealing algorithm which consists of several BFGS 
evaluations. Both the trial designs and the best design obtained are presented. 
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Table 5 Hybrid optimization results 



Reference 

Optimum 

Number of plies 

28 

18 

Root chord (in.) 

15.0 

15.4 

Wall thickness (in) 

0.14 

0.09 

Stacking sequence 
top and bottom walls 

[0°/90°]14 

[07-45°]9 

side walls 

[45°/-45°]14 

[070°]9 

Natural frequency (hz) 

9.4, 34.1,50.7, 116.9, 

8.75,34.7, 46.8, 107.7, 


154.7, 163.6 

159.5, 170.1 

Flutter point 

29 hz, 2nd mode 

74 hz, 5th mode 



Figure 25. Penalty function iteration history 
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ref. opt. 

Figure 26 Weight reduction through hybrid optimization 

Initially, the flutter constraint is violated which results in very large values of the penalty 
function which are not presented due to the scale of the graph. Due to heuristic nature of 
the optimization algorithm, a large number of function evaluations are required which 
make it computationally expensive. The optimum state is reached in less than 100 
iterations and the optimization procedure is terminated after 250 iterations since no better 
design could be found. There is a significant reduction in the weight of the structural 
member of the wing (32 percent, Figure 26) along with a large improvement in the flutter 
speed (75 percent. Figure 27) due to the optimization. The Tsai-Wu stress criterion is 
satisfied by the reference as well as the optimal design (Figure 28). Since the wing root 
chord for the reference and the optimal wings are nearly same (Table 5), weight reduction 
is due to the fewer number of plies in the optimal wing. Through optimization of the 
stacking sequence, even a lower wall thickness provides higher flutter speed. 
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ref. opt. 


Figure 27. Flutter speed constraint for hybrid optimization. 



ref. opt. 


Figure 28. Tsai-Wu failure criterion for hybrid optimization 

Study of the frequencies and modes for flutter show important trends (Table 6). For 
the reference wing, the second mode, which is a coupled bending and lag mode with a 
natural frequency of 34 Hz, flutters at 29 Hz. The first torsion mode is the sixth mode 
with a natural frequency of 164 Hz. At the flutter condition, the frequencies of the 
second and the sixth modes almost coalesce and the sixth mode also flutters at a slightly 
higher speed. For the optimal wing, flutter occurs at the fifth mode (at 74 Hz), which is 
the first torsion mode with the natural frequency of 160 Hz. The first four modes are 
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bending or coupled bending/lag modes. Thus optimization essentially stiffens the 
bending modes to increase the flutter speed. 

Since the aspect ratio and taper ratio are fixed in this study, smaller root chord also 
means smaller span and lower wing area. The optimal value of the wing root chord is 
very near to the minimum value specified to achieve higher stiffness. This trend is 
expected in the absence of other design considerations such as wing loading (which 
affects landing/take-off speed, maneuverability etc.) and internal fuel volume. 

4,4 MDO of Composite Wings 

The developed MDO procedure with refined analysis and optimization technique 
(gradient bassed) is applied to the design of a composite wing. The wing configuration of 
a high speed business jet type airplane is selected for the optimization. Though the 
developed procedure is capable of handling multiple objective functions (and 
constraints), the example considered here includes only one objective function (with two 
constraints). The Laplace domain method of flutter analysis is used in this optimization. 

4.4. 1 Optimization Problem 

The objective is to minimize the weight (W) of the box beam which represents the 
load carrying structural member of the wing. The cruise Mach number and altitude are 
assumed to be 0.85 and 26000 ft, respectively. In accordance to established design 
practice, the flutter/divergence speed is required to exceed 1.20 times the cruise speed. 
Therefore, the flutter/divergence dynamic pressure (qf) is constrained to be at least 3.7 
psi. Constraints are imposed on the maximum allowable stresses due to static loading of 
the wing corresponding to a 3g maneuver. An elliptic spanwise load distribution is 
assumed. The Tsai-Wu failure criterion [84] is imposed on the ply stresses at the root 
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section. The objective function and constraints are expressed in a similar manner as in 
the previous section (Equations (64) - (67)). The composite ply angles and laminate 
thicknesses are used as continuous design variables. 

4.3.2 Results and Discussions 

The selected reference wing has aspect ratio of 10, taper ratio of 0.35 and mid-chord 
sweep of 30° (Figure 29). The wing geometry is held fixed during optimization. The 
width of the box beam is assumed to equal half the wing chord. The wing thickness-to- 
chord ratio is assumed to be 12 percent, which accommodates the box beam with height 
to width ratio of 20 percent. The box beam is constructed from Carbon-PEEK composite 
material [81] whose properties are presented in Table 4. 


Y 



Figure 29. Wing geometry for MDO of composite wings 
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The box beam walls are assumed to be made of eight unidirectional composite 
laminates (each consisting of several plies of identical orientation), stacked 
symmetrically. It is assumed that the laminates in all four walls have the same thickness 
and that the thicknesses do not vary spanwise or chordwise. The top and the bottom 
walls have the same ply angles, as is the case for the two side walls. The ply angles and 
laminate thicknesses are allowed to vary, thus resulting in 12 design variables. The ply 
angles are allowed to vary between -90° and +90° and the laminate thicknesses vary 
between 0.006 inch and 0.048 inch during optimization. A reference configuration with 
conventional ply orientations of (0°/90°/+45°/-45°)s for each wall and uniform laminate 

thickness of 0.024 inch is used as the initial design. 

The structural analysis of the box beam is performed with 10 elements spanwise and a 
single element chordwise. The first six normal modes of vibration are used in the 
subsequent analysis. The mode shapes are normalized such that the generalized masses 
are unity. The generalized aerodynamic forces (GAP) are computed at M=0.85 for 20 
values of the reduced frequency between 0 and 1. The wing is divided into 48 panels. 
Rational function approximation of the GAP's is performed using four denominator 
coefficients. These coefficients are selected to be 0.2, 0.4, 0.6 and 0.8 which produce a 
very good fit. Convergence is reached in 10 cycles during which the objective function 
and the constraints are evaluated 208 times. The total time taken is 2 hr. 50 min. on Sun 
Ultra 1 workstation. Results obtained using the optimization procedure are presented in 
Table 6 and in figure 30. 


Table 6. Optimization results for MDO of composite wings 



Reference 

Optimum 

Stacking Sequence 

top /bottom walls 

(0°/90°/+45°/-45°)s 

(-8°/707287-13°)s 
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side walls 
Laminate thickness 
(all walls, in.) 

Wall thickness (in.) 
Natural frequency 
(first six modes, Hz) 

Wing structural weight (lb) 

Flutter dynamic pressure 

(psi) 

Tsai-Wu failure criteria 


(0°/90°/+457-45°)s 

(0.024/0.024/0.024/0.024)s 

0.192 

1.34,5.55, 15.20, 
22.63,31.94,37.94 

129.2 

3.0 

1.32 


(157487-57-24°)s 
(0.01 9/0.023/0.024/0.023)s 

0.178 

1.58, 6.45, 17.44, 
22.30,35.79,37.51 

119.1 

3.7 

0.89 



Reference Optimum 


4-/I 



Reference Optimum 


(a) 


(b) 
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(c) 


Figure 30. Optimization results for MDO of composite wings 


The structural weight of the wing reduces by eight percent from 129.2 lb. to 1 19.1 lb. 
(Figure 30(a)), whereas the flutter dynamic pressure increases by 23 percent from 3.0 psi 
to 3.7 psi (Figure 30(b)). Also, from Table 6 and Figure 30(c), it is observed that the 
Tsai-Wu stress criterion is strongly violated in the reference design (Fc=1.32). Through 
optimization of the stacking sequence, even a lower wall thickness (resulting in reduced 
weight) satisfies the stress constraint (Fc=0.89). The ply angles for the optimized wing 
are (-8°/70°/28°/-13°)s in the top and bottom walls and (15°/48°/-5°/-24)s° in the side 

walls. It should be noted that the buckling constraits are not imposed in this 
optimization. The laminate thicknesses remain almost unchanged, except for the surface 
laminates in each wall. It is to be noted that in a practical design, standard ply angles and 
thickness are used. Therefore, the nearest values available must be used. 



73 


Study of the mode shapes (Figures 31-34) indicates similarity between reference and 
optimum designs. Modes 1, 2, 3 and 5 are bending modes and modes 4 and 6 are 
pitching modes. There is very little coupling between the two types of modes. However, 
the natural frequencies (Table 6) indicate an important change from reference to optimum 
wing. The frequencies for the bending modes increase by about 15 percent, whereas 
those for the pitching modes decrease slightly. As discussed later, it is the bending 
modes that exhibit particular flutter characteristics in this wing. 



Figure 31 . Bending mode shape for reference composite wing 
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Figure 34. Pitching mode shape for optimized composite wing 

Flutter is often caused when the frequencies of any two modes coalesce and the 
damping of either of the modes goes to zero. For divergence, the frequency also drops to 
zero along with the damping. Figures 35-40 present root loci and frequency and damping 
history for the reference and the optimized wings. In both cases, frequencies of modes 1 
and 2 coalesce and the damping of mode 1 becomes unstable. The optimized wing has 
larger modal damping and somewhat wider separation between these two modes, leading 
to higher flutter dynamic pressure. 
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Figure 35. Dynamic pressure root locus for reference composite wing 



Figure 36. Frequency vs. dynamic pressure for reference composite wing 












4.5. MDO of Smart Composite Wings 


The developed MDO procedure is applied to the design of a smart composite wing, 
with PZT actuators / sensors bonded to top and bottom surfaces of the box beam. The 
wing configuration is similar to a long range high speed business jet. The dimensions are 
scaled down to reduce the number of finite elements and aerodynamic panels. The 
optimization problem now is associated with multiple objective functions and constraints. 
All design variables are treated as continuous variables. 

4.5.1 Optimization Problem 

The objective is to minimize the structural weight of the wing while maximizing the 
aerodynamic efficiency (lift-to-drag ratio). Wing sweep, wall thiekness of the box beam, 
ply orientations and the thickness of PZT actuators / sensors are design variables. The 
problem is formulated as a minimization problem and the two objective functions are as 
follows. 

Minimize: fj = W 

fj^D/L (71) 

where W is the weight of the box beam (including actuators), D is the drag and L is the 
lift. The cruise Mach number and altitude are assumed to be 0.90 and 30000 ft, 
respectively. It is required that the aircraft does not encounter flutter/divergence 
anywhere within the flight envelope. Therefore, a constraint is imposed on the flutter 
dynamic pressure (qf ) which is stated as follows. 
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This constraint implies that the flutter dynamic pressure must not be less than 8.0 psi 
(M=0.9, SL). The wing strength constraint is based on a 3g maneuver at the cruise Mach 
number and altitude. The use of wing loading, suitable for this class of aircraft (70 Ib./sq. 
ft.), and the selected wing area, the ultimate load for design (ultimate load factor = 4.5) is 
determined. The spanwise load distribution is obtained from the panel code [71] 
described earlier. A strength constraint is imposed on the ply stresses at the root section 
where material failure is most likely to occur. The Tsai-Wu failure criterion (discussed 
earlier in section 4.3) is used and this constraint is stated as follows. 

g2=F,-l<0 (73) 

The box beam walls are assumed to be made of eight unidirectional Graphite-Epoxy 
(Table 1) composite laminates (each consisting of several plies of identical orientation), 
stacked symmetrically. In an effort to reduce the computational expense, the following 
assumptions are made about the ply orientations. The top and the bottom walls are 
assumed to have identical ply orientations. The ply angles for the two side walls are also 
maintained identical (but may be different from the top/bottom walls). Therefore, eight 
design variables are used to completely describe the stacking sequence. 

Top/bottom walls: (G] / 02 / 83 / 04)^. = (<|>i <|>4\ (74) 

Sidewalls: (6i/02^®3'^®4)s “ 

where the ply angles are allowed to vary between -90 deg. and +90 deg. It is assumed 
that the laminates thicknesses do not vary spanwise or chordwise and all four walls have 
the same thickness. The laminate thickness design variables are defined as follows. 

(t,/ 17/ 13/14)^ = <l>12)s 
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where the thickness of each laminate can vary between 0.01 in. and 0.03 in. The mid- 
chord sweep of the wing is defined as a design variable with the upper and lower bounds 
specified as follows. 

25»<+,j<40» (77) 

It is assumed that a pair of actuators / sensors are surface bonded to the top and the 
bottom of the box beam (covering the area of the first element along the wing root). The 
thickness of the actuators / sensors (tp), assumed to be identical, is used as a design 

variable. 

tp=(j)i4 0.005" < ^,4 <0.015" (78) 

4.5.2 Results and Discussions 

The reference wing has root chord of 50 in., semi-span of 1 60 in., aspect ratio of 
9.5, taper ratio of 0.35, thickness-to-chord ratio of 1 1 percent and mid-chord sweep of 30 
degree. The wing sweep is varied during the optimization and all other parameters are 
held constant. The length of the box beam and the side wall areas change as wing sweep 
is varied. The width of the box beam is assumed equal to half the wing chord. A 
reference configuration with conventional ply orientations of (0°/90°/+45°/-45°)s for each 

wall and uniform laminate thickness of 0.02 in. (wall thickness = 0.16 in.) is used as the 
initial design. The structural characteristics of the PZT actuators / sensors (G-1 195) [48], 
which affect mass and stiffness of the wing, as well as the induced strain are included in 
the analysis. The applied voltage is such that the maximum electric field for PZT 
materials (15000 volts/in.) is not exceeded. 

The structural analysis of the box beam is performed with 1 0 elements spanwise and a 
single element chordwise. The first six normal modes of vibration are used for the 
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subsequent analysis. The mode shapes are normalized such that the generalized masses 
are unity. The generalized aerodynamic forces (GAP) are computed at M=0.90 for 20 
values of the reduced frequency between 0 and 1 . The wing is divided into 48 panels. 
Rational function approximation of the GAP's is performed using four denominator 
coefficients (equation (49)). These coefficients are selected to be 0.2, 0.4, 0.6 and 0.8 
which produce a very good fit. 


Table 7. Optimization results for MDO of smart composite wing 



Reference 

Optimum 

Stacking Sequence 

top /bottom walls 

(0790°/+457-45°)s 

(-17907327-44°)s 

side walls 

(07907+457-45°)s 

(27907447-36°)s 

Laminate thickness 

(0.02/0.02/0.02/0.02)s 

(0.025/0.01 1/0.01 3/0.012)s 

(all walls, in.) 

Wall thickness (in.) 

0.160 

0.122 

Mid-chord sweep (deg.) 

30 

39 

PZT actuator thickness (in.) 

0.005 

0.011 

Applied voltage (volt) 

75 

165 

Natural frequency 

1.64, 7.56,21.84, 

1.46, 6.53, 18.43, 

(first six modes, Hz) 

33.29, 46.35, 52.37 

27.27, 38.09, 44.02 

Wing structural weight (lb) 

58.62 

47.28 

Lift-to-drag ratio 

11.24 

14.71 

Plotter dynamic pressure 

6.9 

8.0 

(psi) 

Tsai-Wu failure criteria 

0.79 

0.71 
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The developed MDO procedure is computationally efficient as each run takes 

approximately one minute on Sun Ultra 1 workstation for the numerical example 

presented. Convergence is reached in 12 cycles during which the objective function and 

the constraints are evaluated 259 times. Figure 41 shows the reference and the optimum 

wing geometries. The optimization results are summarized in Table 7 and Figure 42. 

The structural weight of the wing (box beam with actuator / sensor) reduces by 1 9 percent 

(from 58.62 lb. to 47.28 lb.) and the lift-to-drag ratio increases by 31 percent (from 1 1.24 

to 14.71). The flutter dynamic pressure for the optimized wing increases from 6.9 psi to 

8.0 psi and, therefore, satisfies the design constraint. The Tsai-Wu strenght constraint is 

well satisfied by both the reference and the optimized designs. 

Y 



Figure 41 . Reference and optimized geometry for smart composite wing 




Figure 42. Optimization results for MDO of smart composite wing 
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The reduction in the weight is achieved through lower wall thickness, which reduces 
from 0.160 in. to 0.122 in. Even with lower thickness of the box beam walls, due to the 
optimized ply orientations, laminate thicknesses and (increased) actuator thickness, the 
flutter dynamic pressure increases to the required value in the optimum configuration. 
These design variables actually contribute more to the flutter dynamic pressure than is 
apparent from the final results. The higher sweep angle of the optimized wing (39 deg.) 
has a detrimental effect on flutter, but it reduces the compressibility drag. Though the 
higher wing sweep produces larger induced drag, the total drag reduces to enhance 
aerodynamic efficiency of the wing substantially. Study of the ply orientations shows 
relatively small changes in the positive angle plies (+45 deg. to +32 deg.) for the 
top/bottom walls and negative angle plies (-45 deg. to -36 deg.) for the side walls. The 
other ply angles remain nearly unchanged. However, the changes in the laminate 
thicknesses are significant. The thickness of the zero degree laminate increases to 0.025 
in., whereas that for other laminates decreases to about 0.012 in. As shown later, flutter 
is caused by the first two bending modes which are stiffened due to increased thickness of 
the zero degree laminates. It is to be noted that in a practical design, standard ply angles 
and thickness are used. Therefore, the nearest values available must be used. Through 
optimization of the ply angles and laminate (ply) thicknesses, even a lower wall thickness 
(resulting in reduced weight) reduces the maximum stress in the laminates. 

Study of the mode shapes (not shown) indicates similarity between reference and 
optimum designs. Modes 1, 2, 3 and 5 represent bending modes and modes 4 and 6 
represent pitching modes. The natural frequencies reduce from reference to optimum 
design, largely due to the increased wing sweep. The eigenvalues of the characteristic 
equation vary with dynamic pressure. A root locus constructed by varying the altitude for 
a given Mach number yields flutter dynamic pressure when one of the real roots becomes 
zero and the imaginary root is non-zero. The variations in frequency and modal damping 
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with dynamic pressure provide an insight into the onset of flutter/divergence. Flutter is 
often caused when the frequencies of any two modes coalesce and the damping of either 
of these modes goes to zero. For divergence, the frequency also drops to zero along with 
the damping. Figures 43-48 present root loci and frequency and damping histories for 
reference and optimized wings. In both cases, frequencies of modes 1 and 2 coalesce and 
the damping of mode 1 becomes unstable. The optimized wing has larger modal 
damping leading to higher flutter dynamic pressure, although the frequency separation is 
smaller between these two modes. 



Figure 43. Root locus for smart composite wing (reference design) 
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Figure 44. Frequency history for smart composite wing (reference design) 
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Figure 45. Damping history for smart composite wing (reference design) 
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Figure 46. Root locus for smart composite wing (optimum design) 
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Figure 47. Frequency history for smart composite wing (optimum design) 



5. MDO of Turbomachinery Blades 


Aerodynamic and heat transfer design objectives are integrated along with various 
mechanical constraints on the blade geometry for the multidisciplinary design 
optimization of turbine blade profile. The K-S function method is used for multiobjective 
optimization. The methods used for blade geometric modeling and aerodynamic and heat 
transfer analyses are briefly described. A numerical example is presented showing the 
benefits of application of the developed MDO procedure. 

5.1 Blade Modeling 

Bezier curves have properties that make them ideal choice to represent complex 
shapes [65]. These curves can be explicitly expressed through the use of Bernstein 
polynomials. The representation of airfoil geometry by Bezier-Bemstein polynomials has 
been successfully used for shape optimization [61, 64]. A two-dimensional boundary is 
defined by Bezier-Bemstein curve of degree n as follows. 

b"(t)= i biBj’Ct) (79) 

The vector of Bezier control points bj, consisting (n+1) values of x and y coordinates of 
the control points, are varied to generate different Bezier curves. The n-th degree 
Bernstein polynomials are given by 

B;(t) = -^t'(l-t)"-J (80) 

■' j!(n-j)! 


where the values of interpolation points t lie between 0 and 1 . 
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This approach can accurately represent a complex shape with relatively small number 
of geometric control points (Figure 49), which are used as design variables. Hence, a 
reduced number of design variables may be adopted in the design optimization procedure. 
The Bezier curves also possess high order of continuity and the end points pass through 
the first and last control points. 



Figure 49. Bezier-Bemstein representation of turbine blade geometry 
5.2 Analysis 

5.2.1 Aerodynamic Analysis 

Since the results of the optimization procedure depend on the accuracy of the analysis 
techniques used, it is important to integrate reliable analyses within the optimization loop. 
The complex flow environment in turbomachinery components can be best described by 
three-dimensional analysis. However, the coupling of such codes within a closed loop 
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optimization procedure can be computationally very expensive. In the current research, 
the aerodynamic analysis is performed using the RVCQ3D (Rotor Viscous Code Quasi-3- 
D) code [66-67]. It is a rapid code capable of analyzing blade-to-blade viscous flows in 
turbomachines. The analysis is based on the thin shear layer approximation of the 
Navier-Stokes equations. The flow equations are mapped to a body-fitted coordinate 
system and a periodic C-shaped grid is used. Second-order finite differences and an 
explicit multistage Runge-Kutta scheme are implemented for time marching flow 
solution. For turbulent flow calculations, the Baldwin-Lomax turbulence model is 
available. The method includes the quasi-three-dimensional effects of rotation, radius 
change and stream surface thickness variation. The code has been validated for several 
test cases and has been used for many applications [85]. 

Grid generation is accomplished using the GRAPE (GRids about Airfoils using 
Poisson's Equation) code [68-69]. This code performs grid generation by solving Poisson 
equations with arbitrarily specified inner and outer boundary points. The desired grid 
spacing and intersection angles at the boundaries are obtained through proper choice of 
the forcing terms in the Poisson equations. 


5.2.2 Heat Transfer Analysis 

The temperature distribution in the blade interior [86] is obtained by solving the 
following equation of two-dimensional heat conduction. 


dx 


. 3T 


- 4 -- 


dy 


, dT 

k — 

. 


= 0 


(81) 


In the above equation, T is the local blade temperature and k is the thermal conductivity 
of the blade material. The finite element method is used to solve the boundary value 
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problem. The computational domain is discretized using linear triangular elements. The 
mesh is regenerated as the blade geometry changes during optimization. The temperature 
obtained from the RVCQ3D solution is used to specify the Dirichlet type nodal boundary 
condition. Using the Galerkin approach, the boundary value problem is reduced to the 
following system of linear equations 

KT = f (82) 

which is solved for the unknown nodal temperatures T . The coefficient matrix K and 
the vector f are evaluated using the finite element formulation. 


5.3 Optimization Problem 

The developed multidisciplinary optimization procedure [87] has been applied to the 
design of a subsonic turbine blade. The objectives are to minimize the total pressure loss 
and maximize the kinetic energy efficiency [85] for the blade-to-blade flow by changing 
the blade profile. The total pressure loss is defined as 
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Since no coolant hole or film cooling is considered in the present design, 
minimization of average and maximum blade temperatures was not included as an 
objective. The blade temperature was treated as a constraint. The maximum and average 
blade temperatures are constrained not to exceed that of the reference blade. Blade chord, 
trailing edge thickness and stagger angle are kept constant during the optimization. 
Accordingly, the two control points at the trailing edge and the leading edge control point 
are not varied. The tangential coordinates of the remaining ten control points are treated 
as design variables, while their axial coordinates remain fixed. The area of the blade 
section is constrained to lie within ten percent of the reference blade area to assure 
structural integrity and to prevent large weight increase. 

5.4 Results and Discussions 

The reference blade (Figure 50) represents a standard section used for turbine design. 
The blade has a finite trailing edge thickness and its chord length and stagger angle are 
0.122 ft and 45 degrees, respectively. The annular cascade is assumed to have 36 blades 
rotating at 2000 rpm. The following flow parameters are specified: Reynolds number = 
6x1 0^ per ft, inlet Mach number (absolute) = 0.21, inlet flow angle (absolute) = 0°, inlet 
total temperature = 1500 K, exit static pressure ratio = 0.70, and relative flow angle at 
trailing edge = -65 degrees. The external flow field around the blade is discretized using 
the GRAPE code with 97 points around the blade and 3 1 points normal to the blade. The 
flow calculations are done with RVCQ3D utilizing four stage Runge-Kutta scheme. The 
blade interior is discretized with approximately 2100 elements for computing the 
temperature distribution. 
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Ref - suction side 



The optimization results are presented in Figures 50-57. The reference and optimized 
blade profiles are shown in Figure 50. Significant changes to the suction and the pressure 
surfaces are observed, though the airfoil area remains unchanged. The optimized blade 
profile represents the most efficient aerodynamic shape at the flow conditions specified. 
In a practical application of the developed MDO procedure, additional constraints may be 
imposed on the blade geometry from other design considerations. The surface Mach 
number and static pressure distributions are shown in Figures 51 and 52, respectively. 
On the pressure side, the small adverse pressure gradient and flow deceleration in the 
leading edge region are eliminated, resulting in a smoother velocity profile. The pressure 
gradient on the suction side is more favorable and hence the flow accelerates to higher 
Mach number near the leading edge. The pressure and the Mach number variations 
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remain non-smooth at the blade trailing edge, because its geometry was held fixed during 
the optimization 
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Figure 52. Static pressure distribution 

The total pressure ratios Pq 2 / Pqi are plotted in Figure 53 for reference and optimized 
blades. The total pressure loss is reduced from 2.9 percent to 1.1 percent. It should be 
noted that even a small reduction in the turbine loss has a significant effect on the overall 
performance of the engine. There is a substantial increase (5.5 percent) in the kinetic 
energy efficiency due to optimization (Figure 54). The temperature distribution in the 
blade interior are presented in Figures 55 and 56 for reference and optimized profiles, 
respectively. The blade temperature constraint is satisfied as shown in Figure 57, where 
the temperature of the optimized blade is normalized with respect to that of the reference 
blade. 



Reference Optimum 


Figure 53. Exit total pressure ratio 





6. Concluding Remarks 


A new multidisciplinary design optimization procedure has been developed for the 
conceptual design of composite wings with surface bonded piezoelectric actuators / 
sensors. The analysis and optimization methods used are computationally efficient and 
sufficiently rigorous for using the developed MDO procedure for actual design 
applications. The optimization procedure for smart composite wing design involves the 
coupling of structural mechanics (including smart material), aeroelasticity and 
aerodynamics. The load carrying member of the wing is idealized and represented as a 
composite box beam. Each wall of the box beam is analyzed as a composite laminate 
using a refined higher-order displacement field to account for the variations in transverse 
shear stresses through the thickness. Detailed structural modeling issues associated with 
piezoelectric actuation of composite structures are included. This structural model is 
suitable for analyzing both thin- or thick-walled constructions. The governing equations 
of motion are solved using the finite element method to analyze practical wing 
geometries. 

The wing steady and unsteady aerodynamic loads are obtained using a panel code 
based on the constant-pressure lifting surface method. This method utilizes the linearized 
aerodynamic potential theory for compressible flows. Two methods for 
flutter/divergence analysis have been implemented. The V-g method has been used to 
predict accurate flutter/divergence speed. The Laplace domain method of flutter 
prediction involves approximating generalized aerodynamic forces, but it produces root- 
loci of the system which give an insight into the physical phenomena leading to 
flutter/ divergence . 
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In the gradient based optimization procedure, the Kreisselmeier-Steinhauser (K-S) 
function approach is used to efficiently integrate the objective functions and constraints 
into a single envelope function. The resulting unconstrained optimization problem is 
solved using the Broyden-Fletcher-Goldherg-Shanno algorithm. For optimization 
problems involving both continuous and discrete design variables, a procedure has been 
developed using the hybrid optimization technique. Based on the several numerical 
examples presented, the following important observations are made. 

1 . Flutter solution is obtained for an example wing using the present approach and 
the classical laminate theory (CLT), which neglects transverse shear stresses. The 
flutter speed using CLT is about 10 percent higher compared to that using the 
present approach. This significant difference in flutter prediction between the two 
theories can be a very critical issue since the flutter phenomenon is catastrophic in 
nature. The difference in the flutter speed is obviosly due to the presence of 
through-the-thickness transverse shear stresses which are ignored by CLT. This 
example establishes the significance of the refined higher-order displacement field 
on the aeroelastic stability of composite wings. 

2. The effect of composite ply orientations on flutter/divergence dynamic pressure 
has been studied, using the Laplace domain method. It has been shown that the 
existence of various coupling modes, for different ply orientations, strongly 
influence the wing aeroelastic characteristics. 

3. In the numerical example for the hybrid optimization technique, wing root chord 
and wall thickness are used as continuous design variables, whereas the ply 
orientations are treated as discrete variables. The wing weight is used as objective 
function which is minimized with constraints on flutter/divergence speed and 
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stresses at the wing root due to the specified air loads. The optimized design 
reduces the wing weight by 32 percent while satisfying the constraints. However, 
due to the use of simulated armealing within the optimization algorithm, a large 
number of function evaluations are required. This makes the procedure 
computationally expensive. 

4. The developed MDO procedure with refined analysis methods and sophisticated 
(gradient bassed) optimization technique has been applied to the design of a 
composite wing. The wing configuration of a high speed business jet type 
airplane is selected for the optimization. The optimization converges in 10 cycles 
taking only 2 hr. 50 min. on Sun Ultra 1 workstation. The optimized design has 
significantly lower wing weight (eight percent) and higher flutter dynamic 
pressure (23 percent). The wing strength constraint, though severely violated in 
the initial design, is met by the optimized design. 

5. The higher-order theory based composite box beam model has been extended to 
include piezoelectric actuators / sensors bonded to top and bottom surfaces. The 
optimization problem is formulated with the objective of simultaneously 
minimizing wing weight and maximizing its aerodynamic efficiency. Design 
variables include composite ply orientations, ply thicknesses, wing sweep, and 
piezoelectric actuator thickness. Constraints are placed on the flutter/divergence 
dynamic pressure and wing root stresses. The maximum electric field applied to 
the actuators is restricted by the coercive field of the PZT material. The 
developed MDO procedure is computationally efficient as each run takes 
approximately one minute on Sun Ultra 1 workstation for the numerical example 
presented. Convergence is reached in 12 cycles during which the objective 
function and the constraints were evaluated 259 times. The optimal design shows 
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physically meaningful results. The weight of the wing is reduced by 19 percent 
and the lift-to-drag ratio increases by 3 1 percent. The flutter dynamic pressure for 
the optimized wing increases from 6.9 psi to 8.0 psi and the Tsai-Wu strenght 
constraint is well satisfied. The wing sweep is increased to 39 degree for reduced 
compressibility drag (hence increased aerodynamic efficiency). The thickness of 
zero degree plies is substantially increased which enhances the damping of 
bending modes leading to higher flutter dynamic pressure. The larger thickness of 
the PZT actuators also helps increase flutter dynamic pressure. 

The development of a new multidisciplinary optimization procedure for the design of 
turbine blades is presented. The procedure integrates aerodynamic and heat transfer 
design objectives along with mechanical constraints on blade geometry. Bezier-Bernstein 
representation of the blade profile leads to a relatively small set of design variables. 
Viscous blade-to-blade flow is calculated through thin shear layer approximation of the 
Navier-Stokes equation. A Poisson’s equation based grid generator provides the grid for 
the flow solution. The maximum and average blade temperatures are obtained through a 
finite element analysis. Total pressure loss is minimized and the exit kinetic energy 
efficiency is maximized with constraints on blade temperatures and geometry. The K-S 
function approach is used to solve the multiobjective constrained nonlinear optimization 
problem. The results for the numerical example show significant improvements after 
optimization. The total pressure loss is reduced by 1.8 percent and there is 5.5 percent 
increase of the kinetic energy efficiency. The maximum and average blade temperatures 
for the optimum blade are lower than the reference case. Other design constraints such as 
airfoil area, chord, trailing edge thickness and stagger angle are satisfied. 



7. References 


1. Venkayya, V. B., "Introduction: Historical Perspective and Future Directions," 
Structural Optimization: Status and Promise, edited by M. P. Kamat, Vol. 150, 
Progress in Astronautics and Aeronautics, AIAA, Washington, D. C., 1993, pp. 
1 - 10 . 

2. Ashley, H., "On Making Things the Best - Aeronautical Uses of Optimization," 
Journal of Aircraft, Vol. 19, No. 1, Jan. 1982, pp. 5-28. 

3. Miura, H. and Neill, D. J., "Applications to Fixed-Wing Aircraft and Spacecraft," 
Structural Optimization: Status and Promise, edited by M. P. Kamat, Vol. 150, 
Progress in Astronautics and Aeronautics, AIAA, Washington, D. C., 1993, pp. 
705-742. 

4. Duysinx, P. and Fleury, C., "Optimization Software: View from Europe," 
Structural Optimization: Status and Promise, edited by M. P. Kamat, Vol. 150, 
Progress in Astronautics and Aeronautics, AIAA, Washington, D. C., 1993, pp. 
807-850. 

5. Johnson, E. H., "Tools for Structural Optimization," Structural Optimization: 
Status and Promise, edited by M. P. Kamat, Vol. 150, Progress in Astronautics 
and Aeronautics, AIAA, Washington, D. C., 1993, pp. 851-864. 

6. Sobieszczanski-Sobieski, J. and Haftka, R. T., "Multidisciplinary Aerospace 
Design Optimization: Survey of Recent Developments," AIAA 96-0711, Proc. 
34th Aerospace Sciences Meeting and Exhibit, Jan. 1996, Reno, Nevada. 

7. Kroo, I., "MDO Applications in Preliminary Design: Status and Directions," 
AIAA 97-1408, Proc., 38th AIAA/ASME/ASCE/AHS/ASC Strctures, Structural 
Dynamics, and Materials Conference, April 7-10, 1997, Kissimmee, Florida. 



105 


8. Mason, W., "Analytic Models for Technology Integration in Aircraft Design," 
AIAA Paper 90-3262, Sept. 1990. 

9. Wakayama, S. and Kroo, I., "Subsonic Wing Planform Design Using 
Multidisciplinary Optimization," Journal of Aircraft, Vol. 32, No. 4, 1995, pp. 
746-753. 

10. Barthelemy, J. -F. M., Wrenn, G. A., Dovi, A. R. and Coen, P. G., "Integrating 
Aerodynamics and Structures in the Minimum Weight Design of a Supersonic 
transport Wing," AIAA 92-2372, Proc., 33rd AIAA/ASME/AHS/ASC Strctures, 
Structural Dynamics, and Materials Conference, April 13-15, 1992, Dallas, Texas. 

11. Grossman, B., Haftka, R. T., Kao, P. -J., Polen, D. M., Rais-Rohani, M. and 
Sobieszczanski-Sobieski, J., "Integrated Aerodynamic- Structural Design of a 
Transport Wing," Journal of Aircraft, Vol. 27, No. 12, 1990, pp. 1050-1056. 

12. Shirk, M. H., Hertz, T. J. and Weisshaar, T. A., “Aeroelastic Tailoring - Theory, 
Practice, and Promise, “ Journal of Aircraft, Vol. 23, January 1986, pp 6-18. 

13. Johnson, E. H., Venkayya, V. B., "Automated Structural Optimization System, 
Volume I - Theoretical Manual," AFWAL-TR-88-3028, 1988, pp 102 and 177- 
179. 

14. Haftka, R. T., Gurdal, Z. and Kamat, M. P., "Elements of Structural 
Optimization," Kluwer Academic Publishers, Boston, 1990. 

15. Chattopadhyay, A. and McCarthy, T., "Multiobjective Design Optimization of 
Helicopter Rotor Blades with Multidisciplinary Couplings," Structural Systems 
and Industrial Applications, Ed. S. Hernandez and C. A. Brebbia, 1991, pp. 451- 
461 

16. Kreisselmeier, A. and Steinhauser, R., "Systemic Control Design by Optimizing a 
Vector Performance Index," Proc., IFAC Symposium on Computer Aided Design 
of Control SystemA, Zurich, Switzerland, 1979, pp. 1 13-117. 



106 


17. Wrenn, G. A., "An Indirect Method for Numerical Optimization Using the 
Kreisselmeier-Steinhauser Function," NASA CR-4220, 1989. 

18. Fadel, G. M., Riley, M. F. and Barthelemy, J. F. M., "Two-Point Exponential 
Approximation Method for Structural Optimization," Structural Optimization, 2, 
1990, pp. 117-124. 

19. Seeley, C. E. and Chattopadhyay, A., “Development of a Hybrid Technique for 
Continuous/Discrete Optimization,” Proc., 6th AIAA/USAF/NASA/ISSMO 
Symposium on Multidisciplinary Analysis and Optimization, Bellevue, WA, Sept. 
4-6, 1996. 

20. Lynch, R. W., Rogers, W. A. and Brayman, W. W., "Aeroelastic Tailoring of 
Advanced Composite Structures for Military Aircraft," U. S. Air Force Flight 
Dynamics Lab. Rept. AFFDL-TR-76-100, April 1977. 

21. Giles, G. L., "Equivalent Plate Analysis of Aircraft Wing Box Structures with 
General Planform Geometry," Journal of Aircraft, Vol. 23, No. 1 1, 1986, pp. 859- 
864. 

22. Giles, G. L., "Further Generali 2 ^tion of an Equivalent Plate Representation for 
Aircraft Structural Analysis," Journal of Aircraft, Vol. 26, No. 1, 1989, pp. 67-74. 

23. Livne, E., Schmit, L. A. and Friedmann, P. P., "Integrated 
Structure/Control/ Aerodynamic Synthesis of Actively Controlled Composite 
Wings," Journal of Aircraft, Vol. 30, No. 3, 1993, pp. 387-394. 

24. Reddy, J. N., "Energy and Variational Methods in Applied Mechanics," Chapters 
2 and 4, John Wiley & Sons, New York, 1984. 

25. Livne, E., "Recent Developments in Equivalent Plate Modeling for Wing Shape 
Optimization," AIAA 93-1647, Proc., 34th Structures, Structural Dynamics, and 
Materials Conference, La Jolla, CA, April 19-22, 1993, pp. 2998-301 1 . 



107 


26. Livne, E., Sels, R. A. and Bhatia, K. G., "Lessons from Application of Equivalent 
Plate Structural Modeling to an HSCT Wing," Journal of Aircraft, Vol. 3 1, No. 4, 
1994, pp. 953-960. 

27. McCarthy, T. R., "A New Higher-Order Composite Theory for Analysis and 
Design of High Speed Tilt-Rotor Blades," NASA CR 196703, October 1996. 

28. McCarthy, T. R. and Chattopadhyay, A., "A Refined Higher-Order Composite 
Box Beam Theory," Composites Part B 28B (1997), pp.523-534. 

29. McCarthy, T. R. and Chattopadhyay, A., "Investigation of Composite Box Beam 
Dynamics Using a Higher-Order Theory," Composite Structures, 41(1998), pp. 
273-284. 

30. Reddy, J. N., "A Simple Higher-Order Theory for Laminated Composite Plates," 
Journal of Applied Mechanics, Vol. 51, December 1984, pp. 745-752. 

3 1 . Bhatia, K. G. and Wertheimer, J., “Aeroelastic Challenges for a High Speed Civil 
Transport,” A1AA-93-1478-CP. 

32. Chen, P. C., Lee, H. W. and Liu, D. D., "Unsteady Subsonic Aerodynamics for 
Bodies and Wings with External Stores Including Wake Effects," Journal of 
Aircraft, Vol. 30, No. 5, 1993, pp. 618-628. 

33. Bisplinghoff, R. L., Ashley, H. and Halfman, R. L., "Aeroelasticity,” Chapters 3 
and 5, Addison-Weslwy Publishing Co., Inc., Cambridge, 1955. 

34. Karpel, M., "Design for Active and Passive Flutter Suppression and Gust 
Alleviation," NASA CR 3482, 1981. 

35. Roger, K. L., "Airplane Math Modeling Methods for Active Control Design," 
Structural Aspects of Active Controls, AGARD-CP-228, August 1977, pp. 4.1- 


4.11. 



108 


36. Abel, 1., "An Analytical Technique for Predicting the Characteristics of a Flexible 
Wing Equipped With an Active Flutter-Suppression System and Comparison 
With Wind-Tunnel Data," NASA TP- 1367, 1979. 

37. Vepa, R., "Finite State Modeling of Aeroelastic Systems," NASA CR-2779, 1977. 

38. Tiffany, S. H. and Adams, W. M., Jr., “Nonlinear Programming Extensions to 
Rational Function Approximation Methods for Unsteady Aerodynamic Forces,” 
NASA TP 2776, July 1988. 

39. von Karman, T., "Turbulence and Skin Friction," Journal of the Aeronautical 
Sciences, Vol. l,No. 1, 1934, pp. 1-20. 

40. Raymer, D. P., "Aircraft Design: A Conceptual Approach," Chapter 12, AlAA, 
Washington, D. C., 1989. 

41. Shevell, R. S., "Fundamentals of Flight," Chapter 12, 2nd ed., Prentice-Hall, 
Englewood Cliffs, NJ, 1989. 

42. Gallman, J. W., Kaul, R. W., Chandrasekharan, R. M. and Hinson, M. L., 
"Optimization of an Advanced Business Jet," Journal of Aircraft, Vol. 34, No. 3, 
May-June 1997, pp. 288-295. 

43. Chopra, 1., "Review of Current Status of Smart Structures and Integrated 
Systems," Proc., SPIE Vol. 2717, SPIE's 3rd Annual International Symposium on 
Smart Structures and Materials, 1996, pp. 20-62. 

44. Weisshaar, T. A., “Aeroservoelastic Control Concepts with Active Materials”, 
Proc., International Mechanical Engineering Congress and Exposition, Chicago, 
IL,Nov. 1994. 

45. Heeg, J., "An Analytical and Experimental Investigation of Flutter Suppresssion 
via Piezoelectric Actuation", Proc., AIAA 33rd Dynamics Specialists Meeting, 
Washington D.C., pp. 237-247, 1992. 



109 


46. Song, O., Librescu, L. and Rogers, C. A., "Application of Adaptive Technology to 
Static Aeroelestic Control of Wing Structures" AIAA Journal, Vol. 30, No. 12, pp. 
2882-2889, December 1992. 

47. Paige, D. A., Scott, R. C. and Weisshaar, T. A., “Active Control of Composite 
Panel Flutter Using Piezoelectric Materials” Proc., SPIE’s 1993 North American 
Conference on Smart Materials and Structures, Albuquerque, NM, Feb. 1-4, 1993 

48. Crawley, E. F. and Lazarus, K. B., "Induced Strain Actuation of Isotropic and 
Anisotropic Plates" AIAA Journal, Vol. 29, No. 5, June 1991, pp. 944-951. 

49. Chattopadhyay, A. and Seeley, C. E., "A Higher Order Theory for Modeling 
Composite Laminates with Induced Strain Actuators" Composites Part B 
(accepted for publication). 

50. Seeley, C. E., "Analysis and Optimization of Smart Composite Structures 
Including Debonding," Ph. D. Thesis, Arizona State University, May 1997. 

51. Lim, K. B. and Gawronski, W., "Control and Dynamic Systems," Vol. 57, 
Academic Press Inc, 1993, Actuator and Sensor Placement for Control of Flexible 
Structures. 

52. Nam, C., Kim, Y. and Weisshaar, T. A., "Optimal Sizing and Placement of Piezo 
Actuators for Active Flutter Supression," SPIE Vol. 2443, SPIE's 2nd Annual 
International Symposium on Smart Structures and Materials, 1995, pp. 40-5 1 . 

53. Chen, G., Bruno, R. J. and Salama, M., "Optimal Placement of Active/Passive 
Members in Truss Structures Using Simulated Annealing" AIAA Journal, Vol. 29, 
No. 8,pp. 1327-1334, 1991. 

54. Rao, S. S., Pan, T. and Venkayya, V. B., "Optimal Placement of Actuators in 
Actively Controlled Structures Using Genetic Algorithms," AIAA Journal, Vol. 
29, No. 6, pp. 942-943, 1991. 



110 


55. Chattopadhyay, A., Seeley, C. E. and Jha, R., "Aeroelastic Tailoring Using 
Piezoelectric Actuation and Hybrid Optimization," Smart Materials and 
Structures, February 1999. 

56. Zhou, S. W. and Rogers, C. A., “Power Flow and Consumption in 
Piezoelectrically Actuatuated Structures” AIAA Journal, Vol. 33, No. 7, July 
1995. 

57. Glassman, A. J., "Turbine Design and Application," Vols. 1 and 2, NASA SP 290, 
1972. 

58. Wilson, D. G. and Korakianitis, T., "The Design of High-Efficiency 
Turbomachinery and Gas Turbines," Chapter 7, Second Edition, Prentice Hall, 
New Jersey, 1998. 

59. Meauze, G., "Overview on Blading Design Methods," in "Blading Design for 
Axial Turbomachines," AGARD Lecture Series 167, 1989. 

60. Vanderplaats, G. 'bi.," Numerical Optimization Techniques for Engineering 
Design: with Applications," McGraw-Hill, New York, 1984. 

61. Burgreen, G. W. and Baysal, O., "Three-Dimensional Aerodynamic Shape 
Optimization of Supersonic Delta Wings," AIAA-94-4271, 1994. 

62. Chattopadhyay, A., Pagaldipti, N. and Chang, K. T., "A Design Optimization 
Procedure for Efficient Turbine Airfoil Design," Computers and Mathematics 
with Applications, Vol. 26, No. 4, 1993, pp. 21 - 31. 

63. Narayan, J. R.,. Chattopadhyay, A., Pagaldipti, N. and Zhang S., "Integrated 
Aerodynamic and Heat Transfer Optimization Procedure for Turbine Blade 
Design," AIAA Paper No. 95-1479, Proc., 35th AIAA/AHS/ASCE/ASME 
Structures, Structural Dynamics and Materials Conference, New Orleans, 
Louisiana, April 1995. 



Ill 


64. God, S., Gofer, J. I. and Singh, H., " Turbine Airfoil Design Optimization," 
ASME Paper 96-GT-158, Proc., International Gas Turbine and Aeroengine 
Congress and Exposition, June 13-16, 1996, Birmingham, UK. 

65. Farin, G., "Curves and Surfaces for Computer Aided Geometric Design: A 
Practical Guide," Chapter 4, Third Edition, Academic Press, Inc., 1993. 

66. Chima, R. V., "Explicit Multigrid Algorithm for Quasi-Three-Dimensional 
Viscous Flows in Turbomachinery," Journal of Propulsion and Power, Vol. 3, 
No. 5, Sept. - Oct. 1987, pp. 397-405. 

67. Chima, R. V., Turkel, E. and Schaffer, S., "Comparison of Three Explicit 
Multigrid Methods for the Euler and Navier-Stokes Equations," NASA TM- 
88878, January 1987. 

68. Sorenson, R. L., "A Computer Program to Generate Two-Dimensional Grids 
About Airfoils and Other Shapes by Use of Poisson's Equations," NASA TM- 
81198, 1980. 

69. Steger, J. L. and Sorenson, R. L., "Automatic Mesh Point Clustering Near A 
Boundary in Grid Generation with Elliptic Partial Differential Equations," Journal 
of Computational Physics, Vol. 33, No. 3, Dec. 1979, pp. 405-410. 

70. Reddy, J. N., "Mechanics of Laminated Composite Plates, Theory and Analysis," 
Chapter 5, CRC Press, Inc., Boca Raton, Florida, 1997. 

71 . “Documentation of ZONA61 Code, Version C,” ZONA 91-8.3, Zona Technology 
Inc., Mesa, Arizona, October 1994. 

72. Hoadley, S. T., Adams, W. M., Jr. and Silva, W. A., "ISACv5 - Interaction 
Structures, Aerodynamics, and Controls," Version 5.3 User's Manual, NASA TM 
- 100666, April 1994. 

73. Szu, H. and Hartley, R., "Fast Simulated Annealing," Physical Letters A, Vol. 
122, No. 3, June 1987, pp. 157-162. 



112 


74. Chattopadhyay, A., Zhang, S. and Jha, R., "Structural and Aeroelastic Analysis of 
Composite Wing Box Sections Using Higher-Order Laminate Theory," AlAA 96- 
1567, Proc., 37th AlAA/ASME/ASCE/AHS/ASC Structures, Structural 
Dynamics and Materials Conference, Salt Lake City, UT, April 15-19, 1996. 

75. Jha, R. and Chattopadhyay, A., "Development of a Comprehensive Aeroelastic 
Analysis Procedure for Composite Wings Using Laplace Domain Methodology," 
AIAA 97-1026, Proc., 38th AIAA/ASME/ASCE/AHS/ASC Structures, Structural 
Dynamics and Materials Conference, Kissimmee, FL, April 7-10, 1997. 

76. Jha, R. and Chattopadhyay, A., "Aeroelastic Stability of Composite Wings Using 
Higher-Order Laminate Theory," Mathematical Problems in Engineering 
(submitted). 

77. Chattopadhyay, A., Jha, R. and Seeley, C. E. "Application of Hybrid Optimization 
Technique for Improved Aeroelastic Performance of Composite Wings," AIAA 
96-4015, Proc., Sixth AIAA/NASA/USAF Symposium on Multidisciplinary 
Analysis and Optimization, Bellevue, WA, September 4-6, 1996. 

78. Chattopadhyay, A., Seeley, C. E. and Jha, R., "Aeroelastic Tailoring Using 
Piezoelectric Actuation and Hybrid Optimization," Smart Materials and 
Structures, February 1999. 

79. Jha, R. and Chattopadhyay, A., "Multidisciplinary Optimization of Composite 
Wings Using Refined Structural and Aeroelastic Analysis Methodologies," 
Engineering Optimization (in press). 

80. Jha, R. and Chattopadhyay, A., "Smart Composite Wing Design for Optimal 
Aeroelastic Control," AIAA 99-1514, Proc., 40th AIAA/ASME/ASCE/AHS/ASC 
Structures, Structural Dynamics and Materials Conference, St. Louis, MO, April 
12-15, 1999. 



113 


81. Agarwal, B.D. and Broutman, L. J., ''‘‘Analysis and Performance of Fiber 
Composites," John Wiley and Sons, Inc., 1990, p. 437. 

82. E. H. Dowell, (ed.), "A Modern Course in Aeroelasticityf Chapter 3, Kluwer 
Academic Publishers, The Netherlands, 1995. 

83. Weisshaar, T. A. and Foist, B. L., “Vibration and Flutter of Advanced Composite 
Lifting Surfaces,” Proc, AIAA/ASME 24th Structures, Structural Dynamics, and 
Materials Conference, May 1983, pp. 498-508. 

84. Tsai, S. W. and Wu, E. M., 1971, “A General Theory of Strength for Anisotropic 
Materials,” Journal of Composite Materials, 5, pp. 58-80. 

85. Chima, R. V. and Yokota, J. W., "Numerical Analysis of Three-Dimensional 
Viscous Internal Flows," AIAA Journal, Vol. 28, No. 5, May 1990, pp. 798-806. 

86. Talya, S. S., Rajadas, J. N. and Chattopadhyay, A., "Multidisciplinary Design 
Optimization of Film-Cooled Gas Turbine Blades," Seventh AIAA/ASME Joint 
Thermophysics Conference, Albuquerque, NM, July 10-14, 1998. 

87. Jha, R., Chattopadhyay, A. and Rajadas, J. N., "Shape Optimization of 
Turbomachinery Airfoil ," Engineering Optimization (submitted). 



